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I. Introduction 

Many mechanical systems can be modeled in state space format as the 
constant coefficient linear matrix differential equation: 

x * Ax + Bu (1) 

y * Cx + Du (2) 

where x s n- vector of state variables 

u * m-vector of control variables 
y * p- vector of output variables. 

In a typical mechanical system, one would be able to measure many, but perhaps 
not all, of the state variables x. An estimate x of the complete state might 
be obtained through a Luenberger observer or Kalman filter. The control 
system design task then consists of finding an algebraic or dynamic control 
law u(x, y, t) which yields the control signal based on measurable quantities. 

This study is concerned with evaluation of alternative computational 
procedures for obtaining the feedback control law. It is desired to find 
computational methods which 

1) involve only a small number of free parameters (i.e. two or 
three) to be specified by the designer so that minimal user 
interaction or "cut and try" iteration is required, and 

Z) yield robust control , i.e. the controller is insensitive to 

small changes in the A and B matrices and performs satisfactorily 
when the mechanical system is operating away from the nominal 
design point. 

The methods evaluated in this study assume that the full s»,ate is 
measurable, and find a constant m x n feedback matrix K with 

u = Kx. (3) 
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The three methods evaluated are: 

1 

1} the standard linear quadratic regulator (LQR) design method , where 
one minimizes the performance index 

J * j j (x^*Qx + u^Ru)dt 

with Q positive semi-definite and R positive definite matrices of 
appropriate dimension. 

2) minimization of the norm of the feedback matrix, |jKj| via nonlinear 
programming subject to the constraint that the closed-loop eigen- 
values be in a specified domain in the complex plane. 

3) maximize the angles between the closed-loop eigenvectors (or, 
equivalently, make corresponding left and right eigenvectors as 
nearly col inear as possible) in combination with minimizing jjKjj 
also via nonlinear programming subject to the closed-loop eigen- 
value constraint in 2. 

These three methods are called the LQR, min K and robust controller design 
methods, respectively. 

The specific function minimization technique used for design methods 2) 

and 3) is a modification of Powell's conjugate gradient method requiring no 

2 

gradient information. Admittedly, this is not the current state-of-the-art 
in nonlinear programing, but the method is simple and reliable, and adequate 
for this preliminary study. 

The domain of the complex plane chosen for closed-loop eigenvalue place- 
ment in methods 2 and 3 is illustrated in Figure 1. The domain is bounded 
on the right by the maximum real part of eigenvalues, REMAX, and on the left 
by the minimum real part of eigenvalues, REMIN. The top and bottom boundaries 
are specified by the maximum ratio of imaginary and real parts of complex 
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eigenvalues, RTOMAX. The closed-loop eigenvalues are those of the matrix 

A * A + BK (4) 

specified as 

\(A) * {X-j, Xg, »«•» X^/ 

Thus the eigenvalue domain can be specified mathematically as 

REMIN ^ Re (X.) £ REMAX 
1 Im (A 1 )‘ 


Re (x,) 


£ RTOMAX 


(5) 

( 6 ) 


for all i * 1 , 2, . . . , n. 

This choice of eigenvalue domain in the complex plane is based on the 
spectral structure of linear systems^. That is, given any arbitrary control 
function u(t), the time response of system (1) is 

x ( t ) * e At x(o) + ( e A(t>x) Bu(t)dT. (7) 

Jo 

The matrix exponential can be decomposed as 

n \.t u 
e At = £ e 1 v<w,” 

i=l 1 1 

where v^ and are the right- and left-eigenvectors of the A matrix, i.e. 


A = X^v^, Wi H A = X iWi H . 

The eigenvalue problem (8) can be written in matrix form as 

A = MJQ 

where J = diag (x^, \ 2 , x n ) 

M = [v-j v 2 ... v n ] 


( 8 ) 

(9) 


Q = M 


-1 


H 


w. 


H 


In general the eigenvalues will be distinct and stable for the closed-loop 
system. 
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Assume that r of the eigenvalues are real and stable, i.e. 

* -Bj < 0, 1 ■ 1* 2* ...» r 

and that 2c eigenvalues are complex and stable, so that there are c complex 
conjugate pairs 

V+j ■ - ^ + 1lS j- - a j < 0 


x r+c+j “ V+j a j 


l iJj , J — 1, 2, »•>, C' 


Then the matrix exponential can be written as 



where E^ , and G i are the appropriate real projection matrices formed from 
the dyad product of right- and left-eigenvectors. 

As seen from (7) and (10;, the rate of decay of the response x(t) of the 
closed-loop system is characterized by the time constant 


x 


max 





( 11 ) 


A larger time constant x means slower system response, and t is kept from 


becoming large by the right-hand eigenvalue boundary REflAX, i.e. 


Re (x t ) < REMAX implies x < • (12) 
Note that REMAX will always be negative for a stable design. The choice of 
REMAX is governed by (12) to achieve a desired or specified closed-loop time 
constant t. 

The effect of the eigenvalue ratio 


Im (.Vj) 


< RTOMAX 


I Re (Vj) 

on transient response is well known for the standard damped harmonic oscillator 
equation 

x + (2^) x + x = 0. 
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That Is, as shown in Figure 1, RTOMAX and are related by the angle 

* - TAN" 1 (RTOMAX) cos^U). 

In order to have well damped non-oscillatory motion in each mode correspond- 
ing to a complex closed-loop eigenvalue, one can choose ? > .71 or y <45° 
which implies RTOMAX <1. 

The left-hand boundary REMIN of the eigenvalue domain is added only to 
form a closed domain. In general, sending closed-loop eigenvalues far to 
the left in the complex plane requires large entries in the feedback matrix K, 
which is prevented by minimizing | jKj j . However, there is no penalty in terms 
of system stability or transient response if closed-loop eigenvalues have 
large negative real part. 

As stated above, the third or "robust" design method was chosen to yield 
a closed-loop system whose eigenvalues are insensitive to small changes in the 
A and B matrices. The relationship between orthogonality of closed-loop 
eigenvectors and the sensitivity of closed-loop eigenvalues is described in 
the next chapter. 
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II. Eigenvalue Robustness Through Eigenvector Orthogonality 

As previously described (9), the closed-loop system is assumed to 
have distinct eigenvalues C\-j > . x n ) where 

A + BK * A • MJM" 1 
J * diag (X-j, ...,X n ) 
and corresponding (right) eigenvectors 

M * [v^Vg. . * v n 3* 

We assume that the closed-loop system matrix A is perturbed as 

A f E (12) 

due to small changes in A and B, and assess the effects of the perturbations 

5 

E on the eigenvalues {Xj, ...» x n > using first-order perturbation theory. 

Let X and v represent a particular eigenvalue/eigenvector pair (possibly 
complex) with the eigenvector normalized so that 

||v|| - (v H v) 1/2 = 1 

We will proceed to find approximations to an eigenvalue X' and eigenvector 
v' of the perturbed system 

(A + E) v * = XV (13) 

that are near \ and v. Since E is small, i.e. 

1 1 Ej j = 0(e), 0 < e « 1 

the differences \'-x and v'-v will also be small, i.e. 

X' - a * u, j u 1 = 0(e)<< 1 (14) 

V' - V = q, j j q 1 1 = 0(e)«l. (15) 

H 

If v' is also normalized as v' v' = 1, then q will be orthogonal to v, i.e. 

q H v = 0. 

In order to obtain expressions for u and q, let U be any n x (n-1) dimen- 
sional matrix such that [vll] is n x n and unitary, i.e. 

[vU] H [vU] = I. 
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so v is orthogonal to each column of U. The perturbation vector q can then be 
written as a linear combination of the columns of U, i.e. 


Up. 


The final results 


I*' - M < ||E|| ||w H || + ||Ei| 2 1 1 U ( A I - U H AU) _1 U H i 
|X' - X| < e 1 1 W H j | + e 2 1 1 (XI - U H AU f ] \\ 


(16) 


and 


07 ) 


V' - v = (XT - U^U)" 1 U H E v 
| V 1 - V | < e I j (XI - U H AU) -1 j I 
are derived in Appendix E . As equation (16) indicates, for small perturba- 
tions E the length of the left eigenvector corresponding to eigenvalue 
X.j provides a measure of the sensitivity of x^ to variations in A + BK. Since 
the left eigenvectors provide a reciprocal basis for the right eigenvectors, 
it follows that 

U 

= 1 . 


llw-iV 

The angle between vectors w.. anc v^ measured by 


= cos 




( w ^ w .) 1 ' 2 ( v ^ v ,)'/ 2 


(18) 


also provides a measure of the length of w- and the sensitivity of X^ to 
perturbations E. Small angles will indicate that is small, and there- 
fore the eigenvalue X^ is "robust". 

The second result, equation (17), indicates that the sensitivity of 
eigenvector v^ is proportional to the distance between eigenvalue x- and the 

- U- 

rest *if the eigenvalue spectrum of A. Note that U AU has the same eigenvalue 


as A less x^- , i.e. 


\{U H AU> = X{A> 


V 


In order for eigenvector v. to be insensitive to variations in A, we want 
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the eigenvalue spectrum \{A> to be well separated in the complex plane, i. 
no two eigenvalues x ^ should be closely spaced, or |x^ -Xj| should be 
maximized for all i f j. 
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III. Design Problems 

The system matrices used for this design study are fourth order with 

two and three control variables, and represent linearized lateral plane 

aircraft dynamics. The state variables are 

"roll rate (rad/ sec) “ 

roll angle (rad) 
x = 

yaw rate (rad/ sec) 
sideslip angle (rad) 


u 


aileron (rad) 
rudder (rad) 

_Yaw control* (rad) 


*for Hall a/c only. 

The first example is taken from Montgomery and Hatch 8 (1969) and represents 
the lateral dynamics of an early version of the Space Shuttle. The system 
and control matrices are: 


A = 


-0.367984 

1.0 

-0.024209 

-0.258819 


0.0 


0.0 


0.0 


0.017835 


-0.032279 26.18750 

0.267949 0.0 

-0.110395 4.46294 

-0.965926 -0.091072 


B = 


-7.67183 

0.0 

1.96959 


2.06549 

0.0 

-2.33843 


o.o 0.0 


The second example models the lateral dynamics of a T-33 trainer and is 

g 

described in Hall (1971). The yaw control is achieved through asymmetric 
deflection of drag petals mounted on wing tip tanks. 
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- 3.18 


0.0 


0.63 


- 10.6 


1.0 

0.0 

0.0 

0.0 

- 0.06 

0.0 

- 0.27 

4.18 

0.022 

0.0644 

- 0.988 

- 0.151 

- 14.4 

1.5 

1.0 


0.0 

0.0 

0.0 


0.0 

- 2.59 

- 0.96 


0.0 

0.037 

0.0 




IV. Programs NEWSOM and LINEAR 

As described in Chapter I, a matrix K is computed to feedback the full 
state to the control variables as 


u = Kx. (19) 

This chapter describes the program NEWSOM and the algorithm used for the 
minimum K and robust controller design methods. 

The basic purpose of the program NEWSOM (listed in Appendix A) is tc 
minimize an unconstrained function of several variables with inequality con- 
straints imposed as penalty terms added to the cost function. The multivariable 
function is minimized using the Powell conjugate gradient nonlinear programming 
method (Zangwill, 1967) implemented in subroutine POWELL. A key part of the 
nonlinear program is the line search algorithm, performed by subroutine MINPT. 
The line search is performed with a combination of outward stepping with 
doubling of successive step sizes, inward stepping using the golden section 
search routine, and parabolic curve fit. 

The independent variables over which POWELL searches are the nm elements 
of the K matrix. The cost function is defined as the sum of five scalar 

terms f^ each with a weighting factor wT . 

5 

f(x) ■ < 20 > 

1=1 1 1 


The scalar functions are defined as 


^(x) = 


nm o 

s «? 


1/2 


= ! ! 1/ t ( 


K| 


( 21 ) 


f 9 (x) = E max(0, REMIN - Re(.\.)) 2 
L i=l 1 

n ? 

f,(x) = E max(0, Re(.\. ) - REMAX r 

J i=i 1 

A 

f 4 (x) = E max(0, abs(Im(.\.)/Re(.\.)) - RTOMAX)' 
H i=l ' 1 


\ 
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where {\j} are the eigenvalues of the closed loop system A * A + BK. The 
three terms f 2 (x), fg(x), f^(x) constrain the closed-loop eigenvalues to 
remain in the domain r of the complex plane illustrated in Figure 1. The 
fifth scalar function provides a measure of eigenvector orthogonality. 
Note that the eigenvector angles ^ defined here differ from the angle 
introduced in Chapter II. As ^ goes to 90° for all i / j, ^ goes to 
zero. 


+1J = cos 


•i v i\ 




f s (x) - 



i-1 


E - 9o °) 10 

M 


1/10 


The power of ten on each term and the tenth root on the summation are employed 
to achieve equal penalization of small angles >£.j. 

The flowchart of program NEWSOM is presented in Figure 2. As shown, the 
program contains integer flags to control the execution of multiple cases 
(with a separate namelist block for each case) and the amount of printed out- 
put. The program input variables are defined in comment cards at the beginning 
of the program. 

The key program inputs are the initial values of the independent vari- 
ables (the elements of K) and the weighting terms wTl , ..., wT5. To generate 
a minimum K control law, the program is run with wT5 = 0 and all other weights 
chosen to place closed-loop eigenvalues in or near the region r. To generate 
a robust control law, the program is run with primary weight on the term 
f 5 (x) and very little or no weight on the f-j(x) term. 

The second program, LINEAR, was used to compute the LQR controller gains 
for both aircraft. The complete program listing is presented in Appendix B. 
This program was not written specifically for this research project, but was 


"V 


IS 


Figure 2 Flowchart of program NEWSOM 


read namelist NSOMIN from file 
NEWSOM INPUT 



ISKIP > 1 



if IPRINT > 3, then print out 
namelist NSOMIN 


© 


if x(l) = -1, then randomize the 
elements of the initial K matrix, 
K = Cx(j)] 


subroutine RANDS computes random 
numbers -1 < x(j) < 1, j=l, ...nm 


if IPRINT > 1, then print out 
A and B matrices 


compute cost function with initial 
K and print out cost terms, eigen- 
values and angles between eigen- 
vectors 


function COSTF computes cost 
function. If IPRINT > 5, 
print partial output; if IPRINT > 
8, print complete output 


minimize the function f(x) of the 
nm elements of K. 

f(t) = WT1 * f ] ( x) + WT2 * f 2 (x) 

+ WT3 * f 3 (x) + WT4 * f 4 (x) 

+ WT5 * f 5 (x) 


subroutine POWELL minimizes an 
unconstrained multivariable function, 
call subroutines: 

COSTF 

RANDIN 

MINPT 

COSTD 

MINPAR 


















previously developed by an AOE department graduate student, Mr. Mark Hreha. 
Flowcharts, input variable definitions and descriptions of the algorithms 
were not available for LINEAR, but the code is based on a report by Sandell 
and Athans (1974). 


V. Results of programs NEWSOM and LINEAR 

For each of the two design problems, I.e. for the Hall and Montgomery 
aircraft models listed in Chapter III, three minimum K controllers, three 
robust controllers and three LQR controllers were computed. The three 
different controllers for the minimum K and robust cases we obtained by 
specifying different values of the time constant t (or equivalenlty REMAX * -r~)» 
namely x a 1, .5 and .25 sec. Table 1 presents the input data for these cases. 
The values of the weights wTl, ..., wT5 were chosen to provide a good tradeoff 
between the placement of closed-loop eigenvalues in the r domain (controlled 
by wT2, wT3 and wT4), and minimization of ||K|{ (controlled by wTl) or maximiza- 
tion of robustness (controlled wT5). The LQR controllers were generated with 
the performance index weighting matrices Q and R, also listed in Table 1. 

The set of three LQR controllers for each aircraft were obtained by varying 
the control weighting matrix in the performance index as 

R = pl m 
m 

with r = 100. , 1 . , and .04. 

The resulting output K matrices from programs NEWSOM and LINEAR for 
the Montgomery and Hall aircraft are listed in Tables 2 and 3, respectively. 

The corresponding closed-loop eigenvalues and angles <j>.j between eigenvectors 
are presented in Tables 4 and 5 following Chapter VII. 

Once the eighteen feedback cases were generated, they were evaluated 
by comparing trajectory time histories, and by comparing the sensitivity of 
closed-loop eigenvalues to perturbations in the A and B matrices. These 
evaluations are presented in the next four chapters. 
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Table 1 Input data for programs NEWSOM and LINEAR 



Montgomery 


Table 2 Gain Matrices for Montgomery Aircraft 


MIN K 

t * 1.0 

x = 0.5 

r = 0.25 
ROBUST 

t = 1.0 

x = 0.5 

i = 0.25 
LQR 

m = 100 
= 1 

.- * .04 


K MATRIX 


1.4941162E-1 

-4.1928029E-1 

1.0535228E-1 
2. 0296771 E- 2 

1.6827383E0 

3.6583920E0 

-2.3151433E-1 

-4.7954880E-2 

1.7022794E-1 
-6.01 64034E-1 

2. 51 96051 E-l 
-2.7018290E-2 

2. 3499689E0 
4.0587797E0 

-1.1856604E0 

4.8925018E-1 

1.2836323E0 
-3.1 967741 E-l 

8.5040188E-1 

-1.5637993E1 

1 . 3259439E0 
4.4778833E0 

3.2358761 EO 
-1.5113544E0 


K MATRIX 


6.5132E-1 

-3.85575E-1 

3.2589E-2 

-1.92548E-1 

1.60285E0 

3.46768E0 

-7.393E-1 

-6.0737E-1 

6.4359E-1 

-3.8107E-1 

6.764E-1 

-2.44E-1 

2.07949E0 

3.5064E0 

-7.41 E-l 
-4.32244E-1 

1.0368E0 

-3.599E-I 

1.27027E0 

-3.22E-1 

2.05683E0 

3.81589E0 

-7 . 43469E-1 
-5.2096E-1 


K MATRIX 


3.306E-1 

-1.784E-1 

1 .0918E-1 
-3.913E-2 

-6. 9998E-1 
4.8287E-1 

2.6222E0 

-1.4068E0 

1.2583E0 

-1.7087E-1 

1 . 0046E0 
9.53E-2 

-3.068E-1 

1.0189E0 

5.9393E0 

-1.734E0 

5.17E0 

-5.9375E-2 

4.8719E0 

1.1309E0 

1.06E-1 

5.0654E0 

7.2623E0 

-2.5309E0 


Table 3 Gain Matrices for Hall Aircraft 


MIN K 
7* l.o 

t» 0.5 

0.25 

ROBUST 

* = 1.0 

? = 0.5 

r = 0.25 

LQR 
= 100 

= 1 

= 0.04 


K MATRIX 


-3.97326E-2 
-4.35985E-3 
1 .99877E-2 

1.20207E-2 

-9.44925E-1 

2.55934E0 

2.16118E-1 

-3.66113E-3 

4.75163E0 

-9.48512E-1 

-5.47779E-1 

-5.14107E-1 

2.46900E-2 
-1.39400E-3 
-1 .16600E-2 

2.40550E-1 

3.84000E-2 

-4.38G00E-2 

-3.01650E-1 

2.29900E-2 

4.14500E0 

-1.46550E-1 

1.14765E-2 

-3.11300E-1 

3.56954E-1 

-3.94971E-2 

-1.27874E-2 

1.45003E0 

2.40516E-1 

-8.76188E-1 

-1.29224E0 
2. 69001 EO 
3.61001 EO 

-4.19082E0 

-1.24536E1 

1.41500E1 


K MATRIX 


-3.48900E-2 

1.51500E-1 

-2.633Q0EQ 

1 . 37833E-1 
-9.33000E-1 
2.5350QEO 

2. 17500E-1 
2.64400E-3 
5, 53060EQ 

-9.67330E-1 

-6.96500E-1 

-5.25889E-1 

3. 81 900E-2 
-1.62450E-3 
-1 .25480E-2 

4.1 3300E-1 
3.86788E-2 
-4 . 43500E-2 

-2.75920E-1 

2.53255E-2 

4.15544E0 

-1.22950E-1 

7.23800E-3 

-3.05500E-1 

4.99740E-1 
3.09131 E-l 
-3.34400E0 

1.46229E0 

-9.32129E-1 

-5.79063E0 

1 . 1 2387E-1 
-1.13726E-3 
5.97918E0 

-1 .04167E0 
-2.05667E0 
-1 . 36633E0 


K M/T 

•IX 


4 . 51 300E-2 
5.74G00E-3 
3.79800E-3 

9.16000E-2 
2.22000E-2 
1 . 1 6790E-2 

5.69790E-2 

9.24200E-2 

3.63300E-2 

-7.54000E-2 

-2.08000E-2 

-8.92500E-3 

8.65300E-1 

-6.75770E-2 

7.83300E-3 

9.95860E-1 

-5.85950E-2 

1.62460E-2 

1 . 1 7500E-1 
9.39890E-1 
3. 51480E-1 

-5.53190E-1 

-2.79920E-1 

-l.OOOOOE-l 

4.83410E0 

-4.26400E-1 

2.78700E-2 

4.98000E0 
-5.01 500E-1 
1 . 23700E-2 

4. 18090E-1 

4.74800E0 

1.76390E0 

-7.36140E-1 
--3 . 1 5300E0 
-8.60260E-1 
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VI Proqram INTODE 

This FORTRAN program integrates ordinary differential equations and is 
used to give the time trajectory of multivariable linear systems. It also 
provides information on the eigenvalue stability of the linear system. The 
equations describing the system dynamics and the feedback law should be in 
the form 

x * Ax + Bu 
u = Kx 

where x is the n x 1 state vector 
u is the m x 1 control vector 

A is the n x n system matrix 

B is the n x m control matrix 

K is the m x n feedback matrix. 

The program integrates the differential equations using the Runge- 
Kutta 4th order method. In this method for a given differential equation 

^ = f (x,y) , 

we have y. +1 = + §■ (kj + 2k 2 + 2k 3 + k 4 ) 

where k, = iif \x.-, y 4 ) 

k 

k 2 = hf ix,j + jjr, y i + y “) 

k 3 = hf (:^ + y-j + 

k 4 = hf + h, y i + k z ) 

and h is the step size. 

The results of ohe program INTODE are available as a trajectory table 
and two sets of plots, one using the printer and the other using the 
versatec plotter. The program also computes the eigenvalues, eigenvectors 
and the angles between the eigenvectors for the closed loop system. 
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Figure 3 is a basic flowchart of the program INTODE. Input data to 
the program is given 'n the form of a namelist. A number of integer 
variables have been included in the program to provide flexibility in 
program execution. Multiple cases can be run by appropriately setting the 
index ICONT and including multiple cases in the namelist. The amount of 
printout is controlled by the index IPRINT while the index IPLOT controls 
the generation of plots. A complete listing of the program is given in 
Appendix C. A list of variables forming the namelist and their notations 
is included in the program listing. 

The program was used to generate time trajectories for the lateral 
dynamics of the Hall and Montgomery aircrafts. The state and control vari- 
ables in the model used are 
p - roll rate I 

4 > - roll angle ( 

} state variables 
r - yaw rate I 

6 - sideslip angle I 

ea - aileron deflection I 

5r - rudder deflection I 

control variables 

The yaw control is present only on the Hall aircraft. Chapter III gives the 
A and B matrices for these models. 

In this study we are comparing the performance of three types of con- 
trollers. These are the Min K controller, Robust controller and LQR controlle 
For each type of controller we have considered three cases and so there are 
a total of nine cases per aircraft. As mentioned earlier, the Min K and 
robust feedback matrices were generated using the program NEWSON. The LQR 
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Figure 3 Flowchart for program INTODE 



28 







matrices were generated using the program LINEAR. Tables 3 and 2 give these 
feedback matrices for the Hall and Montgomery aircrafts respectively. 
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VII Results from Program INTODE 

First we discuss the results obtained for the Hall aircraft. The time 
trajectories obtained with the program INTODE are presented in figures 4, 

5 and 6. In these plots the time trajectories of bank angle, yaw angle, 
aileron deflection, rudder deflection and yaw control deflection are given 
for the nine control matrices described above. The plots have the angular 
variables in radians vs time in seconds. Initial conditions for these 
plots have yaw and bank angle equal to 0.1 radians and all the other vari- 
ables equal to zero. Some of the main features of the results are given in 
table 4. The desired time response should have the following features: 

(a) the response should settle to the steady state value in minimum time 

(b) the response should have minimum overshoot 

(c) the response should require minimum control effort. 

Cases 1, 2 and 3 correspond to the min K controller. Comparing these 
three cases, we see that the second case (for x = 0.5, where x = -1/REMAX) 
gives the best response. Time required to settle to the steady state value 
is about same for all the three cases but the overshoots are minimum for 
case 2. Also the control effort, given by J"u T udt, is the least for the 
• = 0.5 case. Clearly the third case gives the worst time response. 

The three cases which give the time trajectories for Robust controller 
are cases 4, 5 and 6. The time response for cases 4 and 5 are virtually 
identical to those for cases 1 and 2 respectively. Here too, for x = 0.5 
(case 5) we get the best response, requiring minimum control effort and 
having least overshoots. 

The remaining three cases are for the LQR controller. We notice that 
case 7 (p = 100, where o is defined on page22) requires the minimum control 
effort but has a settl inq time much larger than for the other two cases. 

The best compromise is offered by case 8 (p = 1). 
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Comparison of cases 2, 5 and 8 provides some information on the relative 
performance of the three types of controllers considered in this study. The 
response for case 8 is significantly superior than the response for cases 
2 and 5, which are nearly identical. Thus the LQR controller (with.' * 1) 
provides the best time response for the Hall aircraft. Another interesting 
observation is that moving the eigenvalues further to the left does not 
always improve the time response. As x (or p) decreases, the eigenvalues 
in general move further to the left. When we go from t = 1 (or p = 100) 
to x = 0.5 (or p = 1) the time response improves but when r (or p) is fur- 
ther reduced to 0.25 (or 0.04) the response deteriorates. 

Figures 7, 8 and 9 give the time trajectories obtained for the Mont- 
gomery aircraft. Each plot corresponds to a specific case and contains the 
time response of bank angle, yaw angle, aileron deflection and rudder 
deflection. Initial conditions were same as for the Hall aircraft, i.e., 
bank and yaw angles are taken as 0.1 radians while all other variables are 
equal to zero. These plots have the variables in radians vs time in seconds. 
Table 5 presents the main features of the results. 

Using the criterion mentioned earlier, we may infer that of the three 
cases of Min K controller case 3 (t = 0.25) gives the best time response, 
but this requires a very large maximum rudder deflection (-1.715 radians). 

As this is not practically possible, this case is discarded. Comparing the 
other two cases, we notice that the second case requires slightly more con- 
trol effort but it has a much smaller settling time. Hence of these three 
cases, the best compromise is offered by case 2 (t = 0.5). 

Comparison of the results for Robust controller shows that cases 5 
(r = 0.5) and 6 (: = 0.25) give nearly identical response. Case 6 requires 
slightly lower maximum control deflections and overall control effort, where 
as case 5 has a slightly smaller settling time. If we compare the overall 
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Table 5 Results from the Integration of Montgomery Aircraft Dynamics 
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state deviation (given by JxJ x dt) we find that case 6 gives a. lower value. 
So we may choose case 6 as the one giving the best response among the three 
cases of Robust controller. 

For the LQR controller, we notice that as decreases the time required 
to reach steady state also decreases but the total control effort and maximum 
control deflections increase. For case 9 (r = 0.04) an unacceptably large 
aileron deflection is required (1.213 radians). Case 7 (p x 100) requires 
a large settling time and has a large overshoot in the bank angle response. 

So case 8 (p * 1) has a better time response than cases 7 and 9. 

Of the three cases 2, 6 and 8, case 8 requires the minimum control 
effort and also has the minimum settling time, hence case 8 gives the best 
time response. So for the Montgomery aircraft, the best response is given 
by LQR controller with r = 1.0. Here too moving the eigenvalues further 
to the left (i.e., making the real part more negative) does not necessarily 
mean a better time response. 

f T 

Figure 10 qives the plots of the control effort, ju'u dt. versus state 

f T " 

error ;x x dt for both the aircraft. As shown, the LQR controller yields 

lowest control effort and state error. This is to be expected since the LQR 
is by definition the one which minimizes these integrals. 
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VIII Program PERTB 

This program was written to compare the robustness of the three types 
of controllers considered here. The robustness of a controller is measured 
in terms of the insensitivity of the closed loop eigenvalues to variations 
or uncertainties in the A and B matrices, the more insensitive the closed 
loop eigenvalues, the more robust the controller. 

The elements of the A and B matrices are perturbed around their specified 
value by the help of random numbers, 
i.e. PA (1,0) = (1 + Rand x P) A (I,J) 

where Rand is a random number between- 1 , and 1 

P is a fraction giving the maximum perturbation and PA is the 
perturbed A matrix. 

Similarly PB (1,0) = (1 + Rand x P) B ( I , J) 

The closed loop eigenvalues of this perturbed system are calculated. This 
constitutes a single sample. The program also calculates and stores REMAX, 

REM IN and RTOMAX , 
where 

REMAX is Max (Real part of eigenvalues for a particular sample) 

REMIN is Min (Real part of eigenvalues for a particular sample) 

RTOMAX is Max (Ratio of imaginary to real parts of eigenvalues for 
a particular sample) 

The percentage change in REMAX, REMIN and RTOMAX from the unperturbed 
values are then calculated. The program repeats this for a large number 
of samples, typically 1000. Statistics on REMAX, REMIN and RTOMAX and their 
percentage changes are calculated and printed. The program calculates the 
maximum, minimum, mean and standard deviation and presents the variation 
in the form of tabular histogram. 
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The variation in the eigenvalues due to perturbations in the matrices 
A and B can be presented very elegantly as a scatter diagram on the complex 
plane. The unperturbed eigenvalues are circled so as to indicate the amount 
of scatter in the eigenvalues due to perturbations. Such a pictorial re- 
presentation provides qualitative information on the robustness of the 
system. It is a very helpful tool in comparing the robustness of different 
controllers. 

Figure 11 gives the flow chart for the program PERTB and complete list- 
ing is given in Appendix D . Data is given as the namelist PERT and 
details of the variables forming the namelist are included in the program 
listing. As in program INTODE indices ICONT, IPRINT and IPLOT provide 
flexibility in program execution. Sample program input files are included 
with the program listings in the appendixes C and D. 
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Figure 11 Flowchart for program PERTB 




* 

c 


^Eigenvectors are not required for Monte Carlo runs 
For statistics and plotting 
REMAX = Max(real part of eigenvalues) 

REMIN = Min(real part of eigenvalues) 

RTOMAX = Max(ratio of imaginary to real part of eigenvalues) 
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IX Results from Program PERTB 

In this section we compare the robustness of the controllers. This 
is done by observing the scatter of the eigenvalues due to perturbations 
in matrices A and B. For each case a thousand samples, each with 10 /j per- 
turbations of the elements of the matrices A and B, were taken. 

Figures 12 to 16 give the scatter diagrams of the eigenvalues. Tables 
and 6 present the statistics obtained on REMAX, REMIN and RTOMAX. In the 
scatter diagrams, the unperturbed eigenvalues have been circled. This 
helps in estimating the amount of scatter due to the perturbations. Second- 
ly, a number of perturbed eigenvalues are real and so are plotted on the 
real axis. As such, it is difficult to get an idea about their distribu- 
tion along the real axis. So for the figures, we attach a small randomly 
generated imaginary part to real eigenvalues obtained after perturbation. 

These eigenvalues now form a narrow band around the real axis. The density 
of the band at any location gives an idea of the number of perturbed eigen- 
values on the real axis at that location. 

We first discuss the results obtained for Hall aircraft. On studying 
the plots we notice that for both Min k and Robust controllers, as r decreases, 
the scatter of the eigenvalues obtained after perturbations first decreases 
slightly and then increases. This decrease is more apparent for Min k con- 
troller. The Robust controller gives near identical scatter for both ? =1.0 
and t = 0.5. For any particular value of r, the Robust controller gives a 
little less scatter than the Min k controller. 

As shown by the scatter diagrams, perturbations on A and B matrices 
produce very little change in the eigenvalues of the LQR controller. So 
the LQR controller is much more robust than the Min k and Robust controllers. 

It should be noted that one unperturbed eigenvalue for = 1.0 and two 
unperturbed eigenvalues for o = 0.Q4 are not shown on their respective 
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Figure 16 
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plots as their real parts are less than -10. 

Figures 17 to 21 contain the scatter diagrams for the nine cases of 
the Montgomery aircraft. For both Min k and Robust controllers, the scatter 
of the perturbed eigenvalues continuously increases as r decreases from 
1.0 to 0.5 to 0.25. For any one of the three values of t, both Min k and 
Robust controllers give nearly the same amour t of scatter. 

The scatter obtained for the LQR controller decreases as * decreases 
from 100 to 1 and then to 0.04. The scatter obtained for the LQR control- 
ler is significantly less than that obtained for Min k and Robust controllers. 
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Figure 17 


MONTGOMERY^ fl/C^CLJSED LOOP EIGENVALUE PERTURBATIONS, P». I (MIN K) TAU* 1.0 

3.000 ___ 

l * 


2*000 


1*000 


fc 0. 


-1*000 J 


- 2.000 


■rw> 

‘V ."V?*? «VW 


*b 




. ii‘3 

i 


-3.000 


- 10.000 


-a. coo 


» 




- 6.000 - 4,000 
flCflL CLPHOni 


; 

- 2.000 


0 , 


MONTGOMERY A/C CLOSED LOOP EIGENVALUE PERTURBATIONS, P». I (MIN K)TAU»0.5 

>0/22/60 2 i59 

3. COO , ______ 


2.000 J 


i.ooa 


w o, 
az 
cr 
x 

o 

a: 

x 


-1.000 


-2. COO 


-3.000 



* i 


-to. 000 


- 8.000 


-o. COO 


-4,000 


2.300 0. 


6£flU XflHQflJ 


54 












Figure 19 
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Figure 20 
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Table 7 Statistics on REMAX, REM IN and RTOMAX for Montgomery Aircraft 
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’ X. Conclusions and Recommendations 

Based on the results of this preliminary study, it is concluded that: 

1) Searching directly on the nm elements of the K matrix using non- 
linear programming with penalty terms to impose the closed-loop 
eigenvalue constraints 

\ (A + BK) = r 

^ is a feasible, but not highly attractive method of control system 

design. Feedback gain matrices obtained by this method do not 
i appear to have any advantages over LQR controllers, at least for 

I these design examples. The LQR controllers seem to be naturally 

I 

i more "obust than those obtained from direct pole placement. 

1 2) The hoped for increase in system robustness through orthogonal izing 

| closed-loop eigenvalues was not verified in this study. As shown 

‘ in Chapter II and Appendix E, there is a proven mathematical 

^ relationship between eigenvector orthogonality and eigenvalue 

i sensitivity, but this relationship was not clearly demonstrated 

I 

| for these design examples. 

I 

In this continuing research effort to find new and better methods for 

► 

I robust control system design, it is recommended to proceed in the following 

f 

[ directions: 

1) Consider improving robustness only for the class of LQR controllers. 
Perhaps this can be best achieved by searching directly on the 
diagonal elements of the Q and R matrices via nonlinear programming 
to maximize orthogonality of closed-loop eigenvectors. 

2) Test the controller design methods on a wider variety of mechanical 
systems to evaluate controller characteristics. Certain classes 
of problems (such as linearized aircraft equations of motion) 
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may have unique system and control matrix eigenstructures which 
do not provide a sufficiently general test of the design methods. 
Other prototype design examples which could be considered include: 

a) wing flutter suppression in high-speed aircraft 

b) throttle control of multivariable turbofan engines 

c) mathematical examples containing random system and control 
matrices to establish general mathematical properties. 
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Appendix A NEWSOM Computer Program 


C FllF NEWSOM FORTRAN 
C 

C JULY 22, 1990 
C 

REAL *9 OSFED 

COMMON/ADA TA/.N » M» At 4 , 4 ) » 8( 4, 3 ) , X ( 12 ) 

COMMON/BOAT A/NVARS * REM I N ,REMAX > RTOMAX, WT l,WT2,WT3,WT4,WT5 

COMMON/ I PR INT/ 1 PR INT 

COMMON/NNOUT/NODEV 

COMMUN/ISEED/ISEED 

NAMFL tST/NSOMIN/N,M, FPS, MAXFN , NLOOPS , NSRCH, NODE V , I PRINT, I CUNT, I 0 UT 
l ,P. EM IN, REMAX , RTOMAX , WT 1 , WT2, WT3 , WT4, WT5 , I SKI P , A , B , X 
NAME L I ST /NSMOUT/N»Mr/ A, B» X 
C 

C INPUTS 

C N=OI MENS ION OF THE STATE 

C M=D I MPNS I ON OF CONTROL 

C FPS=CONVERGENCE TOLERANCF OF THE POWELL ITERATION 

C NLOOPS -MAX NO OF POWELL ITERATION LOOPS 

C NSRCH=MAX NO OF JOINTS ON A LINE SEARCH 

C Nf)OFV = S PEC I F I ES THE OUTPUT DEV ICE 1 6=T ERM INA L , 8=F I LE NEWSOM OUTPUT) 

C I PRI NT=CONTROLS THE AMOUNT OF PRINTOUT. 1=N0R MAL PRINTOUT 

C l C3N T= 1 , 1 0 CONTINUE AFTER THtS CASE 

r =0 » TO STOP AFTER THIS CASE 

C I O'J T= l , T 0 WRITE NAMELIST BLOCK IN OUTPUT 

C 0 , NOT TO DO THIS 

C WT1...WT5 =WF I GHT S USED IN COST FUNCTION 

C RE‘11 N=LEFT BOUNDARY OF E-VALUE OOMAIN 

C R6MAX=RIGHT BOUNDARY OF E-VALUE OOMAIN 

C R TOMA X=MAX I MUM RATIO IM( LAMDA) /RF ( LAMDA ) 

C A=N* N SYSTEM MATRIX 

C B=N*M CONTROL MATRIX 

C X = **N INITIAL FEEn BACK MATRIX 

C 

DATA NCASE/O/ 

CALI >-‘RRSEr (208 ,256,-1 , 1 ) 

I SEED = 566337 

5 RE AO ( 3 , NSOM I N ) 

NCASE = NCASE +- 1 

IK I SKIP.GF. 1 . ANO.ICONT.LE.OJSTOP 
in ISKIP.3E.I ) GO TO 5 
WRITE! NODEV , 100 INCASE 

100 FORM AT ( / / , ' ********** PROGRAM NEwSOM ********** CASE* X 3) 
CALL DAT I ME ( I M , I D, I Y, I H, IMN, AP, * NEWSOM ', NODEV) 

C I SEE:) = I H* 100 + IMN 

I F ( l !>R I NT • GF • 3 ) WRI TF ( NCOEV , NSOM I N ) 

NVARS = N*M 

I Ft X{ l ) .NE.-l. )G0 TO 3 
~ CALL RANDSt ISEED, NVARS, X) 

rtR IT E ( NODFV , 101 ) ISEED, (X( I ) , 1 = 1 , NVARS) 

101 FORMAT (/,* GGURS CALLED WITH ISEED =*,120, 

l /,' RAND NOS =• , 10F7.4,/, { llX, 10F7.4) ) 

DO 6 I = L, NVARS 

6 XII) = 2 • *X ( I ) - l. 

R (^' mi PRINT. GF. DCALL WRTMAT ( A , N , N, 4, ' A •) 


62 


Appendix A NEWSOM Computer Program 


IF(IPRINT.GF.I)CALL WRTMAT(8 f N,M f 4,*B •) 

IPRINT * IPRINT ¥ 8 
/^) FMIN * COSTFIX) 

IPRINT * IPRINT - 8 

(G) CALL POWFLL <NVARS,X. FMIN, EPS, I PRINT, MAXFN.NLOOPS.NSRCH) 
IPRINT * IPRINT + 3 
FMIN * COSTFIX) 

IPRINT = IPRINT - 8 
(Jp IF ( l OUT. GE. I ) WRITE! 7 »NSMOUT) 

IFIICnNT.GF.l )GO TO 8 
STOP 
END 
C 

FUNCTION COSTFIX) 

01 MF NS ION XU) ,01(5,5) ,0215*5) .03(5,6) ,04(5) 

COMPLEX F I G ( 5 ) » 6VEC (5,5) 

COMMON/ AOATA/N,M, A( A, 4) , B( A, 3), XX ( 12) 

CQMMON/BDATA/NVARS,REMIN,REMAX,RTOMAX,WTl,WT2,WT3, WT4,WT5 

COMMON/ I PR I NT/ 1 PRINT 

COMMPN/NCOST/NCOST 

COMMON/NNOUT/NODEV 

COMMON/ El GOAT/ ER( 5),EI( 5) 

NCOS T = NCOST ¥ I 
FL = 0. 

F2 = 0. 

F 3 = 0. 

F4 = 0. 

F5 = 0. 

00 10 1= I , NVARS 
10 FI = FI f X(I)**2 
FI = WTl*SQRT(Fl) 

CALL MULT (B,X,Di,N,M,N,4,M, 3 ) 

CALL MAT ADD ( A , D l , D l , N » N , 4, 5 , 5 ) 

C 00 12 I = I» N 

C DO 12 J = I » N 

C l 2 DM I,J) = X( I)*A(I|J)/X(J) 

IF(I PRI NT. L F . 8 ) GO TO 15 

CALL WRTMAT(X,M,N,M, •FEEDBACK* ) 

CALL WRTMAT(Dl,N.N, 5, • A*9*K •) 

IS IJOR = 1 

CALL ElGRF{Dl»N»5»IJ0B»EI3»EVEC»5»02fIER) 

I F- ( I FR.GT.O ) WRITE (NODEV, 100 ) IE» 

100 FORMAT ( / * FRRCR IN EIGRF IN COSTF IER =',13) 

00 30 1=1, N 

ER( I ) = REAL ( EIG( I ) ) 

FI ( I ) = A I MAG ( E I G { I ) ) 

I C (FR( I ).LT.RFMIN)F? = F? + wT2*( ER ( I J-REMI N) **2 
IF ( ER ( l ) . GT .REMAX )F 3 = F3 WT3* I ER ( I ) -REMAX) **?. 

RATIO = Ai3S ( El ( I ) / ER t I i ) 

IF{RATI0.GT.RT0MAX)F4 = F<+ + WT4*(RATia-RT0MAX)**2 

TF^p = o. 

00 20 J = 1,N 

01( J, I) = R FAL ( E VEC { J , I ) ) 

02(1 , J ) = 0. 

IFfEH t) .LT.O.IOK J»I) = A I MAG ( E VEC ( J, I ) ) 


00 


Appendix A NEWSOM Computer Program 


I F C I.Ea.J)D2{I, I) * 1 . 

20 TEMP * TEMP * DUJ,I)**2 

30 04(1) * SQRT I TEMP ) 

DO 38 I * 2,N 
l 1 * I - L 
00 36 J * l, I l 
TEMP * 0. 

DO 3* K * l »N 

34 TEMP * TEMP + DL (K, 1 }*0L (K* J) 

ARG * TEMP/ ( D4( I)*04(J) ) 

DitI,J) = 57.295 8* ACQS l AM IN l ! A3S ( AR G ) , l • ) ) 

36 F5 * F5 + {90.-D3I I«J) )**10 

IF< I PRINT. GE. 5) 

IWRI TE (NODEV , 104) ( I » J » D3 ( I , J ) , J* l » 1 1 ) 

104 FORMAT ( / , 5( • ANG* , 1 1 * ' » • , 1 1 , ' *«,F7.2I) 

38 CONTINUE 

F 5 * WT5*( F5**. l ) 

I F ( I PRINT. LT.5)G0 TO 60 
I OGT -- 0 

CALL LFQT2F(Ol,N,N,5,D2, IDGT,D3, • 5 

IF (I fr.gt .0)WRITE (NODEV, 109) I ER 

105 FORMAT!/,' FRROR IN LEQT2F IN CP., i IER =',I6) 

DIN * 0. 

02N = 0. 

FFR = 0. 

DO 50 I = 1 , N 
DP 50 J = l.N 
TEMP = 0. 

DO 45 K = 1,N 

45 TEMP = TEMP + D 1 1 I , K) *D2 ( K, J ) 

IF! I . EQ. J ) TFMP = TEMP - 1. 

FRR = ERR + TEMP**2 
DIN = DIN + DU I » J )**2 
50 02 N = D?N + D2 ( I » J ) **2 

FRR = SORT ( frr ) 

DIN = SORT (DIN) 

02N = SQPT(02N) 

60 F = FI + F? + F3 + F4 + F5 
COSTF = F 

IF(IPRINT.LT.5)RETUP.N 
WRITF! NODEV, 106 ) FRR , 0 IN , D2N 

106 FORMAT!/, • IN COSTF ERR, DIN, D2N =',1P3E12.4) 

*R I TE ( NODEV , 1 08 ) REM IN » REMAX » RTQMAX 

103 FORMAT!/, • REMIN - ' ♦ G 10 . 3 , • RFMAX = ',GL0.3,' RTOMAX = • , G L 0 .3) 
WRITE (NODEV, l 07 )WT1,WT2,WT3,WT4,WT5 

107 FORMAT!/,' WTl , wT2, WT3, WT4, WT5 =',6G12.5) 

WRITE(N006V,191 )F1,F?,F3,F4,F5,F 

101 FORMAT ( ' Fl,F2»F3,F4,F5,FT0T =• ,6312.4) 
wRITE! NODEV ,102) (EKJ I) ,I=l,N) 

10? FORMAT!/, • REAL FIG •*' , 5G12.4) 
wR IT E! NODEV , 103 ) ( E H I ) , 1 = 1 , N ) 

103 FORMAT!' I MAG EIG =',5G12.4) 

CPSTF = F 

RETURN 

END 
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********** SUBROUTINE OATIME ********** 

SUBROUTINE DAT I ME ( I MO, I DAY , IYR, I HOURS , I MIN. AMPM, PNAME,NODEV J 
R EAL *8 PNAMF 

DATA AM/ • AM'/. PM/' PM'/ 

CALL DATE! IMO, I DAY, IYR) 

CALL ST IME ( I TIME) 

XHOURS * FLOAT! ITIME)/ lOOOO# 

AMPM * AM 

IF ( XHOURS. GE. 12. )AMPM = PM 
IF ( XHOUR S. GE. 13 . ) XHOURS « XHOURS - 12. 

I HOURS * XHOURS 
XMIN = (XHOURS - I HOURS ) *60 • 

IMIN a XMIN 

IF(N<?DFV.GT.O)WRITE(NODEV, IOOJPNAME, IMO , I DA Y , I YR, I HOURS , I M I N AM PM 
DO FORMAT!/' TIME IN *,A8,' IS 1 » A2, '/ ' , A2, • /• ,A2, 5X , 12, ' : • , 
l I2,3X,A4) 

RETURN 
END 

********** SUBROUTINE RANDS ********** 

SUBROUTINE RANDS ( I X, N, l ) 

DIMENSION Z! 1 ) 

DATA M/1 048576/, FM/ 1040576./, I A/1027/ 

DO 10 I = 1 , N 
I X = MOD ( I A* I X , M ) 

FX = IX 
10 Z ( I ) = FX/FM 
RETURN 
PND 

C FILF POWFLL FORTRAN 

C 

SUBROUTINE PQWEL L ( N , X , F, EPS, (PRINT ,-MAXFNtN LOOPS »NSRCH) 

DIMENSION X ( l ) , DI ST ( 26 ) ,INDX( 25) 

CDMMCN/P0WEL/NN,D(26,26 ) , P ( 26 , 26 ) , NPCOL ,NDCOL 
COMMON/ISEEO/ISEED 
CPMMON/OUTPUT/INFO 
COMMON/NCOST/NCOST 

COMMON/ IMINPT/L IT, ITO, IT I , ITP , IGO, NONU 

C0MMON/NNOUT/NODEV 

DATA DETMIN/.Ol/ 

LIMIT = NLOOPS 
LIMS = NSRCH 
CALL TIMFON 

INFO = I PR I NT 
NCOST = 0 
NONU = 0 
\'N = N 
IT a 1 
Nl = N + 1 
N? = L • 5*N 
F = COSTF ! X ) 

FI = F 
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FOLD * F . 

IF C INFO.GE.l ) WRITE {NODEV, 103 )N, LI MI T ,L IMS , EPS, F , ( X( I) , 1*1 , N) 

103 FORMAT!/' POWELL ITERATION*, 

l /• N, LI MIT, LI MS* ' ,313,*, EPS*',G12.4, 

l /' F *',G14.7, 

l /• X*' ,5F14.»4,/{3X,5F14.4) ) 

IF( INFO.GE.l )WRITE(NODEV, 104) 

104 FORMAT!' IT IS LIT ITO ITI I TP I GO NONU NCOST* ,2X, 'COS 9X, 

l '0IST',10Xf* FRR ' ) 

C INITIALIZE SEARCH DIRECTIONS 

00 l 1*1, N 
P( 1,1) * X! I) 

1 DIST! I ) * 1. 

DlST(Nt) * 1. 

2 NIT * l 
DET * l. 

DO 10 1*1, N 
DO 5 J*i,N 
5 D( J , I ) * 0. 

10 D ( I , I ) * 1. 

GO TO 15 

12 WR ITE ( NOOEV , 105 ) I , DET , DMA X , AL F 

105 FORMAT!* DN1 REJ I , DET, D MAX, ALF*' , 13, 1P3E14.7) 

C SEARCH IN N DIRECTIONS 

15 OMAX = 0. 

C RANDOMIZE THE ORDER OF CHOOSING THE N SEARCH DIRECTIONS 
CALL R AN I NDI N , INOX ) 

IF!IPRINT.GE.3)WRITE!NOOEV, 107 ) I SFFD , ( I NDX( I ) , 1=1 , N) 

107 FORMAT!' I SEED, INDX =* , 1 11 ,2513 ) 

DO 20 11=1, N 
I = INOX! I I ) 

NPCOL = I I 
NDCOL = I 

CALL MINPT!DIST( I ) , F , EPS , L IMS ) 

IF! INFO. GF. 2) 

1 WRITE (NODFV ,100)1, LIT, ITO, ITI, I TP, I GO, NONU,. NCOST, F, DIST! I ) 

100 FORMAT l IX, 15,514, 15, I6.2G13.6) 

IF! I NFG. GF - 2 ) NONU = 0 

IF { A3S ( DIST ( I ) ) .LF. OMAX) GO TO 18 
DMAX = ABS { 0 1 ST ( I ) ) 

IS = I 

18 DO 20 J=1,N 

20 P(J,IIU) = P(J,II) + DIST! I ) *D{ J » I ) 

C FIND NEW SEARCH DIRECTION 

ALF = 0. 

DO 30 J = i » N 

DlJfiNl) = P ( J , N l ) - P(J,l) 

30 ALF = ALF ♦ D!J,N1)**2 

ALF = SQRT(ALF) 

IF! ALF.EQ.O. ) GO TO 90 
C NORMALIZE NFW SEARCH DIRECTION 

DO 40 J=l,N 

40 D ( J , N l ) = 0(J, ND/ALP 

F3 = F2 
c 2 = FI 
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so 


L 0 1 


102 


60 


00 


106 


C 


10 


?0 

100 


Fi * F 
NPCOL * Ml 
NDCOL * N1 
L IMS * 2*NSRCH 

CALL MINPT! DIST(Ni) ,F,EPS,LIMS) 

L I MS * NSRCH 

ERR * AM INI (ABS( Ir~F3)/F) , AQS(F) ) 

00 50 J*l,N 

P t J » 1 ) * P(JtNl) + D I ST { N l ) *D ( J ,N 1 ) 

XIJ) * PU,1) 

IF ( INFO.GG.l) 

l WRITE (NOOEV, 101 ) IT, I $, LIT, ITO, I TI , I TP, IGO»NGNU,NCQS f , F, 01 ST (N 0 , 
1 PRR 

FORMAT! IX, 12, 13,514, 15, l6,3G13.t>) 

IF( I NPO.GF. 1 )NONU * 0 

I F { INFC.GF.4)WRITF(N0DEV, L02) ( XI K ) ,K*l ,N) 

FORMAT ( * X =• ,LP5E13.6,13X,5E13.6) ) 

IFIERR.LE.EPS. AND. IT.GC.2JG0 TO 90 
I F ( IT. GP. LIMIT) GO TO 90 
IFINCOST.GE.MAXFNJGO TO 90 
IT * IT f l 
IF (NIT.0F.N2JG0 TO 2 
NIT * NIT ♦ 1 

r t -iprK Tn <:pf ip tmf .\iPj qcaoi-u nraprripM CNniii n nc nccn 
OETN * OET*DMAX/ALF 

1 T C PETN.LT.DETMIN)GO TO 12 
DO 60 J = 1,N 

01 J, IS) = D ( J , N 1 ) 

D I ST (IS) = DIST(NL) 

9ET = DFTN 
GO TO 15 
CQNT INUE 

CALL T IMECK I NT IME ) 

TIME = NTIME/100. 

WRITE! NOOFV ,106) I T , NCOST, FOLD , F , T I ME , ( X( I ) , l = l » N ) 

FORMAT!/* SEARCH COMPLETE AFTER', 13, ' ITERATIONS AND*, 15, 

1* FUNCTION CALLS. *,/• INITIAL COST = • , IP E 14 . 7 , • , FINAL COST .7, 
l E14.7,' TIME = *,0PF10.5f//' X = ' , 1P5E 14. 7 , / l 4X , 5E 14. 7 ) ) 
RETURN 
END 

SUBROUT INE CO STD ( D I ST , COST , DM IN , CM IN ) 

COMMON/PnwEL/NN, 0(26,26) ,P( 26,26) , NP COL, NDCOL 
CO.M M ON/IMINPT/LIT , ITO, ITU I TP, IGO,NONU 

common/output/info 

COMMON/NNCUT/ NOOEV 
DIMENSION X ( 2 5 ) 

•^0 10 1=1, NN 

XC I ) = P ( I , NPCOL ) + DIST*DU, NDCOL) 

COST = COSTF(X) 

I F {C OST. GF.CM IN) GO TO 20 
CM IN = COST 
DMIN = DIST 

IF{ INFn.GF.4)wRITF(NU0EV, 100 ) I GO, D I ST , COST , DM IN , CM I N 
FORMAT!' 100=*, II,* D 1ST** , 310. 3, • COST=* ,G 14.7 , 
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l « DMTN*’ ,G10» 3, * CMl.N**,G14.7) 

RETURN 
FNO 

SUBROUTINE MINPT(DMIN,CMIN,EPS,L IMIT) 
COMMON/ IMINPT/ IT , ITO, ITI , I TP, IGOrNONU 
COMMON /NNOUT /NOOEV 
COMMON/IPRINT/IPRINT 

DATA OIG»OTAU/i.E10».118O339886/ 

OAT A 0lSTMN/l.E-6/,VBIG/l.E20/ 

D = DMIN 
OZEAQ * 0. 

CMIN = V3IG 
IGO * 3 

CALL COS Tf) ( f.»?. FRO , C » DM I N * CMI N ) 

CD * C 

IF ( ABS (OI.LT.DISTMN)D * OISTMN 
DELTA * A8S(D>/2. 

02 * 0. 

C2 * C 

SUE = S IGN( 1. ,0) 

IT * NO. OF POINTS IN LINE SEARCH 
I TP = NO. OF PARABOLIC FITS 
[TO = NO. OF OUTWARD STEPS 
ITI = NO. OF INWARD STEPS 
IGO * 0 FOR EXIT ON OUTWARD SEARCH 

1 FOR EXIT ON INWARD SEARCH 

2 FOR PXIT ON PARABOLIC SEARCH 

IT * 0 

[TO = o 

ITO = 0 
ITI = 0 

BEGIN THF OUTWARD SEARCH 

IGO = 0 
01 = -310 
03 = BIG 
Cl = VBIG 
C 3 = VBIG 
GO TO 20 

CHANGE THE DIRECTION OF THE OUTWARD SEARCH 

10 I F { D. GT »D2 ) GO TO 12 

01 = 0 
Cl = c 
GO TO 14 
12 03 = D 

C3 = C 

14 IF ( 01. GT.-BIG.AND.03.lt. BIGIGO TO 35 

SIDE = -SIDE 
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C TAKE AN OUTWARD STEP 
C 

20 I TO » I TO + 1 

I p ( I T. GT. I ) DELTA = QELTA*2. 

0 = D2 ' SI D£*DELTA 

I F ( IPRINT.GE.4)WRITE(N0DEV, 100 ) Dl , D2 , D3, C l , C2 ,C3 

100 F0RMAT(/,3X*»DL' ,11X,»D2* ,11X,'D3* ,11X, 'Cl* ,11X,'C2* ,11X,»C3', 

l /, IX, 6G13. 6) 

CALL C0STD(0,C,DMIN,CMIN> 
l F (IT. GE. L IMIT ) RETURN 
IT = IT + 1 
I c ( C. GT.C2 ) GO TO 10 
l F ( D2. GT. D ) GO TO 16 

01 = 02 
Cl = C2 
GO TO 1 R 

16 03 = D2 

C3 = C2 

18 02 = D 

C2 = C 
GO TO 20 
C 

C THIS COMPLETES THE OUTWARD SEARCH, NOW DO INWARD SEARCH 
C 

35 IGO = l 

GO TO 45 

40 I F { MO D ( IT-IPLAST ,3) )45,45, 50 

C 

C DO PARABOLIC FIT 
C 

45 IGO = 2 

DTEST = .7*ABS( D1-D3) 

CALL MINPAR(Dl,02, D3 , C 1 , C2 ,C3 , CO, EPS , ICONV, DMIN ,CMIN) 

I T P = ITP + 1 
IT » IT + 1 
I F ( I T. GE.L I MIT } RETURN 
IF ( ICONV. EQ. 2 ) RETURN 
IF! ABS(Ql--D3) .LE.0TEST)G0 TO 45 
C 

C GIVE UP ON PARABOLIC SEARCH AND GO TO INWARD GOLDEN SECTION SEARCH 
C 

IPLAST - IT 

02 = {D3+DD/2. - SIDE*DTAU*( D3-DI ) 

IGO = 1 

I F ( IP°INT.GE.4)WRITE(N0DEV, l 01 ) D 1 , 03 , C l , C3 

101 FORMAT (/ , 3X , *01 ' ,24X,*D3' , 1 l X , • C 1 ' , 24X , • C 3' , 

1 /»l Y »G13.6,13X,2G13.6,13X,G13.6) 

CALL CO c 5D(D2,C2,DMIN,CMIN) 

I F ( C 2 . GT . AM I N 1 ( C 1 , C 3 ) ) NO NU = NONU + l 
I F ( C2 .GT , AM IN 1 ( C 1 » C3 ) • AND. I PRINT.3E. 4 ) 

1 WRITE! NO DEV, 102)01, D2,D3, Cl, C2»C3 

102 FORMAT (/, ' CAUTION. ..FUNCTION IS NOT UNIMClDAL. 01-3, Cl-3 = 

1 / , IX , 6G10. 3 ) 

r 

C TAKE AN INWARD STEP 
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C 

*50 0 * UmOU/2. 4- SI DE*DTAU*( D3-DI ) 

I F I IPRINT.GE.4)WRITE( NODEV « 100 )DL,D2,D3,C1,C2,C3 
CALL COS TO < 0» Cf OMINfCMIN) 

I F ( C. GT. AMINl (Cl , C3 ) JNONU « NONU 4- l 
I F ( C .GT.AMINUCUC3) .AND. I PRINT. GE. 4) 
l WRITE (NODEV, 103 ) I) l , 02 , D3, D , Cl , C2 , C3 , C 

103 FORMAT ( / * ' CAUTION, . . FUNCTION IS NOT UNIMODAU Dl-3,0, • 
l ,'Cl-3,C *• ,/, IX, AGIO. 3) 

I T I * ITI 4- I 

FRR - AMINU APS(C),ABS( (C-C2J/C2) ) 

IF ( ERR.LE.EPS.AND.C.LT .CQ I RETURN 
IF( IT. GE. LIMIT) RETURN 
IF * IT 4- I 

sine * siGNii., c C2-c ) *$ i de ) 

DT * D 
CT •= C 

IF ( C • GT »C2 ) GO TO 60 
DT = 02 
CT = C.2 
0? = 0 
C2 = C 

60 IFtSIDE.ST.O. )G0 TO 62 

03 = DT 
C3 = CT 
GO TO 40 
62 01 * DT 

Cl = CT 
GO TO 40 
END 
C 

SUBROUTINE MINPAR(Xl»X2»X3,Fl,F2,F3,FU,FP,ICQNV»DMIN»CM[N) 
COMMON/NNOUT/NOOEV 
COMMON/IPRINT/IPRINT 
ICONV = 0 

A l * X 2 - X 3 

42 = X3 - XI 

43 = XL - X? 

DEN = Fl*Al 4- F2*A2 4- P3*A3 
IF ( QEN . E 0.0. ) RETURN 
3 L •- X?**2 - X3**2 
31 - X3**2 - X l * * 2 
33 = X I * * 2 - X2**2 
X4 = .‘>*(Fl*BL4-F2*B?+F3*33)/DEN 
IF{ X4.LF. XI. 0R.X4.GE.X3) RETURN 
I F ( I PR I N T. GF . 4 ) WR I T E ( NOOEV , 100 ) X l , X2 , X3 , F L , F2 « 

100 FORM AT ( / , 3X , ' D l ' , l IX » *02* , 1 1 X, ' 03 * ,1 IX, •Cl' , 1, ? ' * * C? ' » 1 1 X » ' C3 

I / , IX , 6G 1 3. o ) 

CALL COSTO( X 4 , F 4 , DM I N , C M I N ) 

IF(F4.GT.AMINL(Fl,F3) JNONU = NONU 4- l 
IF(F4.GT.AMINI(FI,F3) .AND. I PR INT .GE .4 ) 

1 WRITE(N00EV»IQ1)X1,X2»X3»X4,FL»F2»F3»F4 

101 FORMAT ( / , ' CAUTION. ..FUNCTION IS NOT ON I MODAL • Dl-4,Cl-4 = 

L /, IX, MG 10. 3) 

ER = A MINI (A8S(F4),ABS( (F4-F2 )/F ? ) ) 
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IF (FR.LE*FP«AND.F4.LT . FQ) ICQNV * 2 
IF l P4-F2 ) 10*10*12 
10 XT * X2 

FT * F2 
X2 * X4 
F2 * F4 

IP! X4-XT ) 14 » 14, 16 

14 XI a XT 
F3 = FT 
RETURN 

16 XI “ XT 

FI = FT 
RETURN 

12 IF(X4-X2) 18, 18,20 

IB XI = X4 

FI = F 4 
RETURN 

20 X'3 = X4 

F3 * F4 
RETURN 
END 
C 

SUBROUTINE R ANIND! N, INOX ) 

DIMENSION INDX! l ) , I ND ( 2S ) , RAND! 2S) 

RFAL*8 OSEBD 
COMMON/ I SFFO/ I SFFO 
COMMON/NNOUT/NODEV 
DO in 1 = 1 ,N 

io ino( n = i 

C WRITE (NODFV, 100) I SPED, N,DSEfc!) 

100 FORMAT (/, ' IN RANDINi I S F ED, N, DSEE 0 - ' , I 12 , I 4, IPO 14. 6 ) 
CALL RANDS! I SEED, N, RAND ) 

C WR li P! NODFV *101)N*DSEED* ( RAND! I ) , I=1,N) 

101 FORMAT!/,' AFTER GGUBS IS CALLED, N, DSEED , I *, 1PD14.6, 

l / , • RAND = ' , 1 0 F 8 - 6 ) 

no 20 K=?,N 

I = N + ? - K 
R I = I 

II = I F I X t R I *RANO( I ) ) + l 
I I = MAXO! 1, MINO! I , I I) ) 

INDXI I) * r ND (II) 

IF! I I .FQ.DGO TO 20 
IU = U + l 

do is j=iii,i 

15 INOIJ-1) = IND(J) 

20 CONTINUE 

INDXtl) = INOU) 

RFTURN 

FND 

C F ILF WRTMAT FORTRAN! A l 

C 

C S/S//D 

C FILE OF UTILITY SUBROUTINFS TO SUPPORT VIBE FORTRAN A1 

r 

C WRTMAT - GENERAL MATRIX OUTPUT SUBROUTINE 


\ 
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C 

SUBROUTINE WRTMAT!A,N,M» IA,ANAMF) 

DIMENSION A ( I A» l ) 

COMMON/NNOUT/NPRINT 
RFAL*8 ANAMF 

WR IT E( NPRINT, 100 ) ANAME , N , M 

100 FORMAT!/,* MATRIX » , AS , 3X , • I • , I 3 , • ROwS X*,I3,' COLS)*) 
IFIM.LE. 10)G0 TO 15 

DO 10 I * i , N 

10 WRITE (NPRINT, LOI) C A { I , J) , J=1,M) 

101 FORMAT! /, I IP10E13.5) ) 

RETURN 

15 DO 20 1=1, N 

?0 WRITE! NPRINT, 102) (A! I ,J) , J=1,M5 

102 FORMAT! 1P10F13. 5) 

RETURN 

END 

C 

C TRANSP - TRANSPOSES A MATRIX 

C 

SUBROUTINE TRANSP! A,N1 , N2 *NA) 

DIMENSION A ( N A , 1 ) 

CO.MMON/NNOUT/NOUT 
M = MAX 0 ( N 1 , N2 ) 

IF (M.GT.NA.0R„M„LE.0)G0 TO 90 
Ml = M~ 1 

00 10 1=1, Ml 
11=1+1 
DO 10 J=I1,M 
TEMP = A! I , J) 

A! I, J) = A! J, I) 

ID A! J, I ) = TEMP 

RETURN 

90 WRITFINOUT, 100)Nl, N2,NA 

IDO FORMAT!* *** ERROR IN TRANSP *** N1,N2,NA =*,314) 

RETURN 
END 
C 

C MULT - MATRIX MULT I PL ICA T ION 

C 

SUBROUTINE MULT l A , B , C , N , L , M , N A , NB , NC ) 

01 MENS ION AINA, 1 ) ,B(NB, l) ,CINC, 1 ) 

DOUBLE PRECISION TEMP 

commqn/nnout/nout 

IF(N.GT.MINO(NA,NC) . OR . L . GT . NB ) GO TO 90 
DO 20 I = l , N 
DO 20 J=l,M 
TEMP = 0. 

DO 10 K=l,L 

10 TFMP = TEMP + OBLEIAI I,K) )* 3BLE ! B ( K » J ) ) 

20 Cl I , J ) = TEMP 

RETURN 

90 WRITEINOUT, 100 ) N , NA , NC , L ,NB 

100 FORMAT ( ' *** ERROR IN MULT *** N , NA , NC , L, NB= * , 5 1 5 ) 

RETURN 
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END 

C 

C MSHIFT - TRANSFERSS VECTORS AND MATRICES 

C 

SUBROUTINE MSHIFT(A,B,N*M»NA»NB) 

DIMENSION A!NA,l) ,B(NB, l) 

COMMON/NNQUT/NOUT 
I F ! N.GT.MINO !NA» NB ) ) GO TO 90 
DO 10 I = 1 » N 
DO 10 J= L , M 

10 B( I, J) * A ( I, J) 

RETURN 

90 WRITE! NOUT » 100 ) N , NA , NB 

LOO FORMAT!* *** ERROR IN MSHIFT *** N » NA , NB= ' , 31 5 ) 

RETURN 
END 
C 

C MATADD - MATRIX ADDITION 

C 

SUBROUT INF MAT ADD! A , B , C , N, M, N A , NB , NC ) 

DIMENSION A ( NA » l ) , B ( NB » l ) » C ( NC , 1 ) 

COMMON /NNOUT/NOUT 
IF(N.GT.MINO(NA,NB,NC) )G0 TO 90 
DO 10 1=1, N 
DO 10 J= l , M 

10 C ( I » J ) = A ( ( , J ) + B { I , J ) 

RETURN 

90 WR ITE ( NOUT ,100)N»NA,NB,NC 

100 FORMAT! • *** ERROR IN MATADD ♦*» N ,N A, NB, NC= • , 4 1 5) 

RETURN 
END 
C 

C MAT SUB - MATRIX SUBTRACTION 

C 

SUBROUTINE MAT SUB ( A , B , C, N ,M,N A , NB , NC ) 

DIMENSION A (N A , I ) , 3 ( NB, 1) , C( NC » 1 ) 

COMMON/NNQUT/NOUT 
I F ( N.OT »M INO I NA , NB , NC ) ) GO TO 90 
DU 10 1 = 1, N 
DO 10 J=l,M 

10 C( I , J ) = A ( I , J ) - 3( I, J) 

RETURN 

90 WRITE! NOUT, 100)N,NA,NB,NC 

100 FORMAT!* *** ERROR IN MAT SUB *** N ,NA ,NB, NC=* , 41 5) 

RETURN 
END 
C 

C SWAP - INTERCHANGES TWO VARA&LES 

C 

SUBROUTINE SWAP(A,B) 

C = A 
A = H 
8 = C 
RETURN 
END 


\ 
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C 

C MS MULT - MATRIX*SCALAR MULTIPLICATION 

C 

SUBROUTINE MS MULT (S»A»N»M,NA) 

DIMENSION A ( NAt L ) 

Of 10 I— I f N 
00 10 J=l,M 

10 A ( I , J ) = S*A ( I * J ) 

RETURN 

END 

C 

C ZERO - FILLS A MATRIX WITH ZEROS 

C 

SUBROUTINE ZERO! A, N, M, NA ) 

DIMENSION AINA, 1 ) 

00 10 1=1, N 
DO 10 J=l,M 
10 A( I, J) = 0* 

RETURN 

END 

C 

C 

C 

C IMAT - LOADS AN ARRAY WITH THE IDENTITY MATRIX 

C 

SUBROUTINE IMATI A»N,NA ) 

01 MFNSION A { NA , 1 ) 

DO 10 I = l , N 
00 5 J=l,N 
5 A{ I f J) = 0. 

10 A( I , I) = l. 

RETURN 

FNO 

C 

C 

C ANORM - CALCULATES THE RSS NORM OF A MATRIX 

r 


FUNCTION AN0RMCA,Nl t N2,NA) 

DIMENSION A I NA , 1 ) 

CO M MQN/NNOU T /NOUT 
IF { N 1 • GT . NA )G0 TO 90 
ANORM = 0. 

DO 10 1=1, N1 
00 10 J= 1 , N2 

10 ANORM = ANORM ♦ A(I,J)**2 

ANORM = SORT! ANORM) 

RFTURN 

90 WRITE! NOUT , 100 ) N , NA 

100 FORMAT ( • *** ERROR IN ANORM *** N,NA =*,215) 

RETURN 

C FILE NFWSOM INPUT 

r 

o 

C JULY 22, 1980 

C 

F.NSOMIN 
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N * 4» M * 2, 

EPS * l.E-7, 

MAXFN * 5000 » 

NLOOPS x 40, 

NSRCH * 12, 

NOOFV * 6, 

[PRINT * 1, 

ICONT * l, 

IOUT * 0, 

WT1 * .01, WT 2 “ LOO., WT3 * 500., WT4 * IQO., 

WT5 * 20. , 

REMIN = -?0., RFMAX = -1., PTOMAX = .5, 

A = -0.367984, l., -0.024209, 0.258819, 3*0., 0.017835, -0.03227*) 

0.267949, -0.110395, -0.965926, 26.1S75, 0., 4.46294, -0.091071 
8 = -7. 6718 3, 0., 1.96959, 0., 2.06549, 0., -2.33843, 5*0., 

X =. 1894, -.633, .359, .01327, 2. 366, 4.1414,— 1.313, .6683, 

X = . L494,-. 4193, . 1054, .0203, 1.6827, 

3, 658, -.2315, -.04795, 

X = .6481, -.4288, .1054, -.2476, 1.6827, 

3.638, -.73056, -.04789, 

SEND 
SNSOMIN 
RFMAX = -2. , 

SEND 
CM SOM IN 


REMAX = -4., 
SEND 


SNSDMIN 

M = 3, 

' “3.18, 1«, -0.06, 0.022,3*0., 0.0644, 0.63, 0., -0.27, —0.998 

-10.6, 0., 4.18, -0.151, ' 

8 = -14.4, 3*0., 1.5, 0., -2.59, 0.037, 2*0., -0.96, 0., 

X=-2. 23E- 8 ,-3 • 5E-8, 3. 43E-3 , 1 . 73E- 10 , . 45 , -. 82 , -1 . 8 ,-l .76 , 10. 7 , 73 
X = -3.9 7E— 2 »-4.56F-3, 1.99 F-2, 1.232 E— 2, -.9449, 

2. 5693,. 216 12, -3. 6611 E-3,4. 75 16, -.94851 , 

-.54778, -.5141, 

X = -.03498, .15196, -2.6329, .1373, -.93295, 

2.5353, .2175, 2.6445E-3, 5.5306, -.9673, 

-.6965, -.5259, 

REMAX = -1. , 

SEND 
SNSOMIN 
REMAX = -2. , 

SEND 


,-5.9 


SNSOMIN 
RFMAX = -4., 
ICONT = 0, 
SEND 
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FILE LINEAR FORTRAN 

REAL A( 10,10) ,8(10,10) ,C I 13, 10) ,QI 10,10 ) ,R{ 10,10) 

REAL K ( 10 , 10 ) ,G( l 0,10) , ACL (10,10) ,QII 10,10 ) ,XC( 10 ),UEI 10) 
REAL Al( 10, 10) ,311 10, 10) ,URFF( 10) 

COMMON/ INOU/K IN, KOUT 
COMMGN/MA IN l /NO l M» 0UM1 1 10,10) 

C0MM0N/MAIN2/0UM2I 10,10) 

COMMON/PLOT/ SC ALE, ARRAY (101,11) 

K I N= 3 
K0UT=4 
NO I M= 10 
N=4 
M=3 
I R=4 

CALL MAT 1 01 N , N, A , 4 ) 

CALL EQUATE(N,N,A1,A) 

CALL MATI0IN,M,B,4) 

CALL EQUATE { N,M, 3 1 , 3) 

CALL MAT 10 ( IR,N,C,4) 

CALL MATIOI N, N,Q, 4) 

CALL MATIOI M,M,R,4) 

CALL CDNIN.M, A,8,NCS) 

CALL OBS(N,N,A,C,NCS) 

CALL RES(N»M,A,B*R,Q,K»G»ACL) 

CALL VECT 10 1 N , XC r 4) 

CALL VECT IOI M , UE , 4 ) 

CALL VECTI0(M,UREF,4) 

CALL MATIOI M,N,G, 4) 

CALL LTRACK(AI,N»B1»M»G»XC,UE) 

STOP 
ENO 

SUBROUTINE REG 

SUBROUTINE R EG( N, M , A, B , R , Q, X, G, ACL ) 

01 MENS ION All) ,3(1),R( l) ,QI 1),X( l),G(i) » ACL ( 1) 

DIMENSION RRt 30) ,RII 30) 

C0MMDN/MAINI/N01M»DUMI( 1) 

COMMON /MAIN2/DUM2I 1 ) 

COMMON/ INOU/K IN, KOUT 
IFINOIM.LT. N) WRITEIKOUT, 1000) 

I F{ NOIM.LT.M) WRITEIKOUT, 1000) 

IFINOIM.LT. N) CALL EXIT 
IFINOIM.LT. M) CALL EXIT 
1000 FORMAT!/, 16H DIMENSION ERROR,/) 

CALL ME I GVI N, A ,RR , R I ) 

WRITEIKOUT, 200) 

WRITFIKOUT , 300 ) (RR( I ) , R 1 1 I ) , I*1,N) 

CALL EQUATE I M , M,0UM1 , R ) 

CALL GM INV I M, M , DUMl , DUM2 » MR ♦ 0 ) 

CALL MAT 41 N, m, 0UM2 , B , G ) 

CALL MRIC(N,A,G,Q,X,ACL) 

CALL EQUATE! M ,M ,DUM1 ,R ) 
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CALL GMINV(M, M,0UMl , DUM2 , MR, 3 ) 

CALL MAT6(M,M,N,DUM2, B,DUM1) 

CALL MMUL(DUML,X,M,N,N,3) 

WRITE (KOUT ,60 ) 

CALL MATIOI N»N,X»3) 

WRITE! KOUT ,70) 

CALL maTI0!N*N»ACL»3) 

CALL MEIGV(N* ACL»RR,RI ) 

WRITF(K0UT,400) 

WRITE (KOUT, 300HRRI I ) , RI ( I) » I*l*N) 

WRI TE ( KOUT , BO ) 

CALL MATI 01 M, N» G » 3 ) 

60 FORMAT (// , 19H RICATTI SOLJTION K,/) 

70 FORMAT!//, 31H OPTIMAL CLOSED LOOP MATRIX ACL,/) 

80 FORMAT!//, 22H OPTIMAL GAIN MATRIX G,/) 

200 FORMAT!/, 22H OPEN LOOP EIGENVALUES,/) 

300 FORMAT ( L 3H REAL PART = ,E10.3,13H I MAG PART - ,EI0.3) 
400 FORMAT!/, 27H CLOSED LOOP EIGENVALUES = ,/) 

RFTURN 
END 

SUBROUTINE MRIC 
SURROUTINF MRIC(N,A,S,Q,X,Z) 

01 MENS ION A! L ) , S € L ) ,Q( 1) ,X( 1) , 2 ( 1 ) ,TR!30) ,TIMES( 3 ) 
COMMON/MAINI/NDI M, F( l) 

COMMON/ I NOU/K IN, KOUT 
NO IM1=NDI M+ 1 
TIMES! 1)=.5 
TIMES(2)=2.0 
TIMES! 3)=4,0 
NN=N*NDIM 

T l =- . 5 * ALCG ( XNORM ( N , Q ) ♦ . 000 1 ) 

IF(TI.LT.-L.0)T1=-1.0/T1 
IF(ARS(TL ).LT.l.O)TI=L.O 
T 2= l. *N/ ( L. +XNORM! N, A ) ) *Tl 
T I =T2 
KEY=0 

5 KEY=KEY+I 
10 DO 15 I = l » N 

DO 15 J= I , NN, NO I M 
15 X(J)=-S(J) 

CALL l NT E G ( N » A » X , Z » -T 1 ) 

CALL CACTOR(N,Z,F,MR) 

C POSSIBLE UNCONTROLLABILITY IF MR.NE.N 

IF(MR.EQ.-L) CALL MAT 10 ( N, N , Z, 3 ) 

I F ( MR . EQ.~ 1 ) CALL EXIT 
CALL CM I N V ( N » N , F , Z » MR , 3) 

CALL MAT2(N,N,Z,Z,X) 

C A+SX IS STABLE 

TOL= l . E-5 
ADV= TOL* 1 . E-7 
NN=N*Nni M 
N M l=N- 1 
OH 19 I = 1 , N 
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TR( I 

19 CONTINUE 
TOLl*TOL/lO. 

MAXI T= 30+N 

90 40 IT* I» MAXItT 
CALL MMUHSr X,N»N,N, F) 

CALL MMUL(X ? F,N,N»N,Z) 

00 20 I*L,NN,NOIM 
I I-l+NMI 
00 20 J=t ,11 
X! J)=Ai J)-F! J) 

20 Zt J)*Z(J)+Q(J) 

CALL ML I NEQ ( N , X, Z , X, TUL 1 ) 

L=0 

Cl=9.0 
11 = 1 

00 25 1 = 1, N 

1 FI ASSC X (II )-TR ( I ) ) .LT. ( AQV+TOL*Xt II))) L=L + l 
TR ( I ) = X ( 1 1 ) 

I 1 = 1 IfNDIML 
25 Cl*Cl«-TR! I ) 

l F ( ABS(Cl) .GT.1.F20) GO TO 50 
IFCL.ME.N) GO TO 40 
CALL GMINVINtN, Z,F,.MR,0) 

CALL MMULIS,X,N,N,N,7) 

WRITE! KOUT , 27 ) IT 

27 FORMAT! L7H0RICCATI SOLN IN ,I2,1LH ITERATIONS) 

00 30 I=l,NN,NDlM 

1 1=1 +NMI 

00 30 J= I , I I 
30 Z(J)=A(J)-Z< Jl 

IF(MR.NE.N)WRITE(KOUT ,35)MR 
35 FORMAT (26H0R ICCATI SOLN IS PSD — RANKI 3 ) 

R F TURN 
40 CONTINUE 

WRI TE ( KOUT , 45 ) MAX IT , T 1 

45 FORMAT (26H0R ICCATI NON-CON VE°GENT INI2,LLH I TERAT I ONS , I2H ItfrrrAt- 
I T=F 10. 5 ) 

GO TO 60 

50 WRITE (KOUT » 55) IT»Tl 

55 FORMAT (29H0RICCATI 3L0W JP AT I TER ATI ON I 2 , 1 2H INITIAL T=F10.5) 

60 IF! KEY.EQ.4)CALL EXIT 
T l = T 2*T I MFS ! KEY ) 

WRITF(K0UT,65)TL 

65 FORMAT! 14H0RESET WITH T=F10.5) 

GO TO 5 
END 

SUBROUTINE ML I NEQ 

SUBROUT INF ML INEQI N» A, C, X» TOL ) 

SOLVES A'X+XA+C=0 

A AND X CAN 9F IN SAME LOCATION IF DESIRED 
ANSWFR RETURNED IN C AND X 
0 I ME NS ION A (L ) ,C( l ) , X! 1) ,RR( 30) ,Rtt 30) 
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COMMON/MAIN I /NO I M,F{ l) 
COMMON/MAI N2/Y ( l ) 

COMMON/ 1 NOU/K IN ,KQUT 
NDIMl=NDIM+i 
OT=. 5 
OTi*0 . 

NN»N*NDIM 

00 5 1 1= 1 ,NN , NO I Ml 
5 0T1=0T1+ABS( A{ II) ) 

OT 1=QT 1/N 

IF(DTl.GT.4.0) DT=0T*4.0/DT 1 

1 1 = 1 

no 20 1 = 1 , N 

00 15 JJ=I »NN ,NDIM 
15 Y(JJ)=DT*A{ JJ) 

Y C 1 1 )=YC 1 1 ) — . 5 
20 II=II+NDIM1 

CALL GM INV ( N» N , Y , F » MR » 1 ) 

CALL FQUATF{N,N,Y,A) 

T F ( MR. BQ.N ) GO TO 21 

1 T=0 

DO 18 1= l , NN, NOI Mi 
18 CU) = 1.625 
GO TO 95 

21 CALL MMUL ( C » F , N,N » N , X ) 

C INITIALIZATION OF X,F 

1 = 1 

00 AO 11=1, HN,NDI M 
J=I I 

IFU.EQ.1I GO TO 30 
on 25 JJ=I , II.NDIM 
C ( J)=C( JJ) 

25 J=J+1 
30 ID=J 

00 35 J J=I I , NN, NO I M 
C(J)=DT*OOT(N»F( I I) ,XIJJ) ) 

35 J = J + 1 

F( ID)=F(,I OI + l.O 
40 1=1+1 

50 A0V=T0L*l.F-7 

00 90 IT= 1 » 30 
NFZ = 0 

S I ZE=0 . 0 

CALL MMUL(C,F,N,N,N,X) 

1 = 1 

I 1 = 1 
J=1 

GO TO 70 
60 J=II 

DO 6 5 JJ«l , I I ,NDIM 
C ( J)sC( JJ) 

65 J=J+1 
70 I D= J 

DT1=C( J) 

00 75 J J = 1 1 , NN , NO I M 
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C(J)*C(J)+DOT(N,F(ir ,X(JJ) ) 

75 J*J+l 
J*J~ l 

no 80 JJ*II,J 
80 X ( J J )*F I J J ) 

I F I ABS I C ( ID)-DTI).LT.(ADV+T0L*ABS<CIID)) ) ) NEZ=*NEZ*L 
1*1 + 1 

II*I IfNDIM 

SI ZE = SI ZE+DTl 

IFII.LE.N) GO TO 60 

IF(NEZ.FQ.N) GO TO L 50 

IF (ABS( SIZE) .ST»i*F25) GO TO 95 

CALL MMUL < X, X ,N»N , N» F ) 

90 CONTINUE 
95 WRITEIKOUT » 100) IT 

100 FORMAT ( 33H0L I N EON ALGORITHM NON-CONVEPGENT , l 3 , 10HI TERAT IONS) 
WRITEIKOUT, 1 LO ) 

110 FORMAT I 35HOA .MATRIX EIGENVALUES* .REAL IMAJ) 

CALL MEIGVIN, Y,RR,RI ) 

WRITEIKOUT, 120) (RRI I ) , R 1 1 I) , I* 1 , N ) 

120 FORMAT! 18X,1P2EI2. 3) 

150 CONTINUE 

CALL EQUATE I N , N , X , C ) 

RETURN 

END 

SUBROUTINE INTEG 

SUBROUTINE INTEG I N, A , C , S , T ) 

S= INTEGRAL EA*C*EA • FROM 0 TO T 
C IS DESTROYED 

DIMENSION All), CIL), St I) ,COEF( L5) 

CO M MON/MAIN I/NDIM, XI 1 ) 

NO I M l =^JQ I M+ 1 
NN=N*NDIM 
NMl=N-l 
NT =13 
NTMI=NT- 1 
I N D= 0 

ANORM=XNORMIN, A) 

DT = T 

5 IF IAN0RM*ABS(DT) .LE.0.5) GO TO 10 
DT=DT/2 • 

I ND= I ND+ 1 
GO TO 5 

10 DO 15 I=1,NN,NDIM 
J=I+NM1 
DO 15 JJ=I , J 
15 SI JJ ) = DT*CI JJ) 

TL=0T**2/2. 

DO 25 IT = 3, 17 

CALL MMUL ( A , C » N » N , N » X ) 

DO 20 1=1, N 
I 1 = 1 I- 1 ) *ND I M 
DO 20 J J = I , NN , NO I M 
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II-II+I 

C( JJ)*(X( JJJ+XI 1 1) )*Tl 
20 S( JJ )*S C JJ)*C( J J ) 

25 T1*DT/FL0AT( IT) 

IFIIND.EQ.O) GD TO 100 
C0£F(NT>*1.0 

00 30 I*l,NTMI 
I I*N T-I 

30 COFF(II )*OT*COEF(II+i)/FLOAT( I) 
DO 40 I»1,NN*N0IM 

1 I-I 
J-I+NMI 

DO 35 JJ*I,J 
35 X{ JJ)*A( JJ) *COEF(l) 

X( 1 1 )*X( 1 1 ) +COFF ( 2 ) 

40 1 1 *1 1+NDI Ml 
DO 55 L = 3 » NT 
CALL MMUL(A,X,NtN,N,C) 

I i*i 

T1=C0EF(L) 

OP 55 1= 1 »NN»NDI M 
J= [ +NM1 
00 50 JJ*I,J 
50 X ( JJ )=C f J J ) 

X( II )aX(IIH-TL 
55 1 1 = 1 I+NDIMl 
C X=EXP(A*OT) 

L-0 

60 L = L + i 

CALL MMU L ( X f S » N t N » N » C ) 
r i = i 

00 90 1*1, N 
J= 1 1 

IF ( I.EQ.L) GO TO 75 
00 70 JJ= I t I I ,NQIM 
S ( J J ) =5 ( J ) 

70 J = J+l 
75 00 35 J J= I , N 
K K = J J 

00 30 K= I ,NN, NO I M 
S( J) = S( J H-C(K)*X(KK) 

*0 KK=KK+-NDI M 
35 J= J+ND I M 

OP 97 JJ=I , NN»NDI M 
37 C(JJ)=X(JJ) 

90 ii=i r+NorM 

IFIL.EO.IND) GO TO 100 
CALL MMUL(C,C,N,N,N,X) 

GO TO 60 
100 CONTINUE 
RETURN 
END 

SUBROUTINE GMINV 
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SUBROUTINE GMINV(NR,NC,A,U,MR,MT) 
DIMENSION A(L)tU(L)*S(30) 

common/maini/ndim 

COMMON/INOU/KINf KQUT 

NDI MI*NDI M+I 

TOL* l . E~ 14 

AOV» 1.6-24 

MP *NC 

NRMi*NP.-l 

TOL 1*0 . 

JJ*1 

DO 10 J=i,NC 

S( J)*DOT(NR, A(JJ) ,A( JJ) ) 

IF(SIJ).GT.T0L1)T0L1=S(J) 

10 J J=J J+NDIM 
TOLI=ADV*TOL1 
AOV*TOLL 
J J* 1 

DO 100 J - 1 , NC 
FAC=S( J) 

4Mi=J-l 
JRMsJJ+NRMI 
JCM=JJf JM1 
DO 20 I'-JJ»JCM 
20 U< I)=0. 

U( JCM)»1.0 
IFIJ.EQ. 1) GO TO 54 
KK = 1 

DO 30 K=l,JMi 

IF(SIK) .EQ.l.O) GO TO 30 

TEMP=-DOT INR, A( J J ) ,AIKK) ) 

CALL VADD(K»TEMP»U(JJ) »U(KK) ) 

30 KK=KK+ND I M 
DO 50 L* 1,2 
KK= 1 

DO 50 K=l f JMl 
IF(S(K).EQ.O.) GO TO 50 
TFMP=-OOT (NR* A(JJ) T A ( KK ) ) 

CALL VADD(NRfTEMP,A(JJ ) , AIKK) ) 
CALL VADD(K t TEMP,U( JJ ) ,U(KK) ) 

50 KK=KK + ND IM 

T0L1 = TCL*FAC*-A0V 
FAC=DOT( NR,A( JJ) , A t J J ) ) 

54 IFIFAC.GT.TGLL) GO TO 70 
DO 55 I = J J » J R M 

55 All )=0. 

S( J)=0. 

KK=L 

DO 6 5 <= 1 » J M 1 
I F ( S ( K ) . E'3.0. ) GO TO 65 
TEMP=-DOT ( K , J( KK) ,U( J J) ) 

CALL VAOO<NR,TEMP, A( JJ) ,A(KK) ) 

65 KK^KK+ND I M 

C AC=DGT( J,*Jl JJ) ,U( JJ) ) 

MR=MR-1 
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GO TO 75 
70 S(J)=1.0 
K K = l 

00 72 K-1»JM1 

1 F ( S ( K ) . EQ. 1 • ) GO TO 72 
TEMP=-DOT(NR,A( JJ) ,A(KK) ) 

CALL VADD(K»TEMP» U( JJ ),U(KK) J 

72 KK=KK+ND IM 
75 FAC=1./SQRT(FAC) 

00 80 I=JJ,JRM 
80 A( I) =A( I )*FAC 
DO 85 I=JJ,JCM 
85 U( I)=U( 1)*FAC 
100 J J=J J+ND IM 

IF (MR -FQ.NR.OR.MP.EQ.NC) GO TO 120 
I F ( MT.NE.O ) WR I TE( KOUT* 1 lOJNRf NC#MR 
110 FORMAT ( 1 3» 1HX, I2t 8H M: RANK, 12) 

120 NEND=NC*NDIM 
J J-l 

00 135 J= l , NC 
DO *25 1=1 , NR 
I I = I-J 
SI I ) =0 . 

00 125 KK=JJ,NEND,NDIM 
125 S( I) = S(I)+A{I I«-KK5*UCKK) 

1 I = J 

00 130 1=1, NR 
UIII)=3([) 

130 1 1 = 1 I+NDIM 
135 JJsJJ+NDIMl 
RETURN 
END 

SUBROUTINE FACTOR 

SUBROUTINE FACTOR ( N, A , S , MR ) 

C A=S‘S 

DIMENSION A ( 1 ) , S ( 1 ) 
C0MM0N/MAIN1/NDIM 
COMMON/ I NO U/K IN, KOUT 
ND I M 1=ND I M+ 1 
Mf<=0 

NN=N*NDIM 
T 0L= 1 . E-7 
TOL 1=0 . 

00 l 1=1 , NN,NDI Ml 
R= A3S ( A ( I ) ) 

1 IF(R.GT.T0L1)TQL1=R 
TOL 1=TCL l *1 • E- 12 

1 1 = 1 

00 50 1=1, N 
I M 1=1-1 

00 5 JJ=I ,NN,NOIM 
5 S{JJ)=0. 

1 0= IT -*■ I M 1 
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R=A( TO J-DOT ( I Ml » S ( II) , SI II) ) 

IFIABS(R)«LT.( TOL*A( I 0 ) +TOL I ) ) GO TO 50 
IF (R) 15,50,20 
15 MR*- L 

WRITFIK OUT, 1000) 

1000 FORMAT (37H0TRIED TO FACTOR AN INDEFINITE MATRIX) 
RETURN 

20 SI 10 )=SQRT( R) 

MR=MR+1 

IF { 1 • EQ-. N ) RETURN 
L = I I +NDI M 

00 25 J J=L » NN, NDI M 
I J = J J* IM1 

25 SI I J ) — C A ( I J J-DOT (IM1,S(II),S(JJ)) ) /S 1 1 0) 

50 II=II+NOIM 
RETURN 
END 

SUBROUTINE MEIGV 

SUBROUTINE ME IGVI N, A,RR , R I ) 

DIMENSION A(L),RR!1),RI 1 1 ) , C I 3 1 ) , TEMP 1 30 ) 
CO.MMON/MA INl/NDIM,Xt 1 ) 

NO I.M1=ND I M+ 1 
nn=n*nqim 

DO 1 1=1, N 
DO 1 J=I ,NN,NDIM 
l XI.J)=A(J) 

CIN+l)*l.O 

1 = 1 

5 C1*0,0 

DO 10 1= 1 ,NN, NDIM1 
10 Cl=Cl-X( I ) 

C l=C 1/FLOAT ( L ) 

I = N+ 1-L 
C I I ) - C 1 

DO 15 1=1 ,NN»N0IM1 
15 X ( I)=X( I ) +C 1 

IF(l.FQ.N5 GO TO 50 
00 40 I = 1 , N 
J J = i 

00 35 J= I ,NN ,NDI M 

r i=o. 

XK= J- I 

DO 2 5 K= I ,NN, .NDI M 
KK'=K K+ I 

25 Cl=CUX(K)*AIKK) 

TEMPI JJ )=C1 
35 JJ=JJ+l 
J J=1 

DO 40 J= I ,NN, NDIM 
X(J)=TEMP(JJ) 

40 JJ=JJ+l 
L = L+ I 
GO TO 5 
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50 CALL P0LRT(C,N,RR,RI) 

RETURN 

END 

SUBROUTINE POLRT 

SUBROUTINE POLRT I A, N, U, V ) 

DETERMINES THE ROOTS OF AN N-TH ORDER POYNOMI AL 
X**N * A( N)*X**(N-l } ♦ ... + A ( 1 ) * 0 
WHERE: U(N) WILL CONTAIN THE ROOT REAL PARTS 

V(N) WILL CONTAIN ROOT IMAGINARY PARTS 
AIN) WILL BE DESTROYED OURING COMPUTATION 
DIMENSION Ml).U(l)tVll) 

1 NR = N 

10 I F* I NR-2 ) 61,71,11 

11 IF (AID) 20,12,20 

12 UINR ) = 0.0 

V ( NR) = 0.0 
NR = NR - 1 

CALL VECTFQINR, A( 2) ,A( 1) ) 

GO TO 10 
20 EMIN = l. 

TOL •= .1 
V4 = 1.0 
P = 0. 

Q - 0. 

R = 0.0 

30 UINR) - A (NR) - P 

UINR-1) = AINR-l) - P*U ( NR ) - Q 

V I NR ) * UINR) - P 

VINR-U = U(NR-l) - P+VINR) - Q 
I = NR-2 

35 UII) = AID - (P*U(I + 1) + Q*U(I+2)) 

VII) = UII) - (P*V(I+l) + Q*V( 1+2 ) ) 

I = I-l 

IF I I.GT.O) GO TO 35 

40 IF (AI2)) 42, +1,42 

41 F = U I 2 ) / A ( 1 ) 

GO TO 43 

42 F = U ( 2 ) / A I 2 ) 

43 F = AMAXl I APS I E ) , 1.0E-6)*AMAX1(ABS(U( 1)/A{ 1) ) , l.QE-6) 

IF ( E.LE. l.E-12) GO TO 70 

IF (E. GE. EMIN) GO TO 44 

THIS FORCES EMIN TO HOLD STFAOY FOR 5 ITERATIONS 
EMIN = E 
TOL = FMIN*0.7 
GO TO 45 

THIS WILL ALLOW AN ERROR X*EMIN ONLY AFTER N ITERATIONS 
WHERE X = I 1 • l ) **N 

44 IF ( E.LT.TOL ) GO TO 73 

45 CBAR = VI 2 ) - U I 2 ) 

IF (NR.GT.3) V4 = VI 4) 

0 = V(3)**2 - C3AR*V4 
IF (0) 47,46,4? 

46 P = P - 2.0 
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0 = q*(Q+i.o) 

GO TO 50 

47 P = P +■ ( ).)( 2 ) *V ( 3 ) - U(l)*V4)/D 

Q * 0 «■ (-Um*CBAR * U(1)*V(3))/D 
50 UINR) = A (NR' + R 
V I NR ) = U( NR ) + R 

1 = NR - l 

55 U( I ) = A C I ) R*U( I* L * 

VII) = IM I ) + R*V ( I -t* 1 ) 

I = 1-1 

IF I I.GT.O) GO TO 55 
F = ABSCUI i) /All)) 

IF IF.LE.l.E— 12) GO T(J 60 
IF (6.GE.EMIN) GO TO 56 
EMIN = E 
TOL = EMI N*0. 7 
GO TO 57 


56 

IF 

IF. 

LT 

.TOL) 

GO TO 60 

57 

IF 

1 V 1 2 ) 

. NE. 0 ) 

GO TO 58 


R = 

R+l. 




GO 

TO 

59 



59 

R = 

R 

- 

U( t) /VI 2) 

59 

TOL 

= 

TOL*L. L 



GO 

TO 

30 




C STORF A SINGLE RFAL ROOT 

60 CALL VECTEQ(NR-1,U(2) ,A) 

GO TO 62 

61 R = -AIL) 

62 UINR) = R 

V I NR ) = 0.0 
NR = NR- 1 
GO TO 80 

C STORE A PAIR OF ROOTS 

70 CALL VECTEQ I NR-2 , U( 3 ) * A ) 



GO 

TO 72 


71 

P = 

A( 2) 



Q = 

All) 


72 

P = 

(-0.5) *P 


0 = 

P*P - 

Q 


IP 

10) 75, 

78 

75 

UINR) = P 



UINR-l) = P 
V ( NR ) = -SQRTI-D) 
VINR-l) = — V I NR ) 

GO TO 79 
73 VI NR) = 0*0 
VINR-l) = 0.0 
0 = ABSIP) + SORT I 0 ) 
IF (P.LT.O.O) 0 = -0 
UINR) = D 
UINR-1) = Q/0 

79 NR NR- 2 

80 IF INR.GT.O) GO TO 10 
RETURN 
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SUBROUTINE FQUATE 

SUBROUTINE EQUATE ( NR, NC , A , B ) 

DIMENSION All), BID 
COMMON /HA I N 1 /NO I M 
NN=NC*NDI M 
NR1=NR-1 

DO 1 J=1,NN,NDI.M 
1 1 =J +NR l 
DO 1 IJ=J, I I 
A ( I J )=B ( IJ) 

1 CONTINUE 
RETURN 
END 

SUBROUTINE M AT4 

SUBROUTINE MAT4IN 1 ,N2 ,X, Y, Z ) 

Z-YXY* X=X' IS N2XN2, Y IS N1XN2, l IS N1XN1 
DIMENSION X(l) ,Y(1) ,Z( 1) 

COMMON/MAI N 1 /NO IM 
CALL MMUL(Y,X|NI,N2»N2»Z) 

NN2=N2*NDIM 
DO 5 1=1, Ml 

I M 1= I- 1 

I I = 1 M 1 *ND I M 
J J=I +1 l 
DO 3 J= I , N 1 
TEMP=0. 

KK=J 

DO 1 K= I ,NN2 , NO I M 
TE>M P = TEMP+Y ( K ) *Z ( KK ) 

L XK=KK+MDIM 
Z( JJ) =TEMP 
3 J J = J J*ND I M 
JJ = [ 

K=I 1+1 
KK= I I + IM1 
DO 5 J=K,KK 
Z ( JJ )=Z ( J ) 

JJ=J J+NDIM 
5 CONTINUE 
RETURN 
END 

SUBROUTINE MAT6 

SUBROUTINE MAT 6 ( N1 , N2 , N3 , X , Y , Z ) 

DIMENSION X( l ) ,Y( l) ,Z( 1) 

C COMPUTE Z = XY ' WHEN X IS N1XN2, Y IS N3XN2, Z IS NLXN3 

COMMON/MAIN1/NDIM 
DO 2 1=1 , N 1 
DO 2 J = 1 , N 3 
TM=G. 
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00 l K=1,N2 

1 TM=TM+X( H-(K-n*NDIM)*Y(J + (K-l)*NDIrt) 

2 Z(I* ( J-i)*NDIM)=TM 
RETURN 

END 

SUBROUTINE MMUL 

SUBROUTINE MMUL ( X, Y, N l , N2 , N 3, Z ) 
DIMENSION X ( 1 ) ,Y( I ) ,Z(l) 

C0MM0N/MA I N 1/NDI M 
NEND3-ND IM*N3 
NEND2=N0IM*N2 
00 I 1*1, N1 

00 I J=I ,NEND3, NDIM 
TM=0. 

K* l 

K K= J- I 
5 KK=KK+1 

TM=TM+X ( K ) * Y( KK ) 

K=K+NQ I M 

I F ( K . L E. NEND2 ) GO TO 5 
I Z(J)=TM 
RETURN 
END 

SUBROUTINE MATIO 

SUBROUTINE MAT I 0( NR , NC» X, 10) 

DIMENSION X(I) 

COMMON/MAINL/NOIM 
COMMON/ INOU/K IN,K0UT 
JEND=NC*NDIM 
IE(I0.EQ«4) GO TO 40 
I p ( I 0 . EQ. 3 ) GO TO 20 
5 CONTINUE 

1 NPUT 

DO 10 1=1, NR 

READ (KIN, 1000 H X( I J ) , I J= I , JEND , NDI M) 

10 CONTINUE 

IFflO.EQ.DRETURN 

OUTPUT 

20 DO 30 1=1, NR 

WRITE! K0U1 , 1001) ( X( I J) , I J= I , JEND, NO I M ) 
30 CONTINUE 
GO TO 50 

TITLE 

40 READ (KIN, 1002 ) 

WR ITE( KOUT , 1003) 
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WRITEIKOUT »1002) 

WR ITE! KOUT ,1004) 

GO TO 5 
50 CONTINUE 

1000 FORMAT!8EiO.O) 

1001 FORMAT! IX* IP10E13.4) 

1002 FORMAT! IX, 79H 

C ) 

1003 FORMAT!//) 

1004 FORMAT!/) 

RETURN 
END 

FUNCTION DOT ( NR» A» B) 

DIMENSION A! 1) ,B(1) 

D0T=0. 

DO 1 1=1, NR 
1 DOT=UOT+A{ I )*B< I ) 

RETURN 
END 

SUBROUTINE VADD 

SUBROUTINE VADD I N , C 1 , A , B ) 

DIMFNSION A! I) ,3(1) 

00 i 1=1, N 
1 A! I )=A ( I )+Cl*Bl I ) 

RETURN ■ 

END 

FUNCTION XNORM! N, A ) 

COMPUTES AN APPROXIMATION TO NORM OF A — NOT A BOUND 
DIMENSION All) 

COMMON/MA I Nl/ND I M 
ND I M 1=ND I M+ 1 
NN=N ♦NO I M 
CL=0. 

TR=A ( 1 ) 

I F ( N« EQ« 1 ) GO TO 20 

1 = 2 

DO 10 I I = ND I Ml , NN »NDIM 
J= 1 1 

DO 5 JJ=I , I I , ND IM 
C1=C l+ABS i A ! J ) *A ( J J ) ) 

5 J=J + 1 

TP=TR+AI J) 

10 1 = 1+1 

TR = TR/FLQAT f N ) 

DO 15 I I = 1,NN,NDIM1 
15 C l=C 1+ ( A ( II)~TR)**2 
20 XNOR M=A8S(TR)*SQRT!Cl) 

RETURN 
END 

SUBROUTINE MAT2 


SUBROUTINE m A T2 ( N 1 , N2 , X , Y , Z ) 
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Z = XY» X,Y=Nl*N2,Z=Z» 

Z AMD Y CAN 3E EQUIVALENT 
DIMENSION X( L) ,Y(1),Z(1) 

COM MON /MAI N 1 /NDI M 
NDIMl=NDIM+l 
NN2=N2*NDIM 
11=1 

00 10 1 = 1 , N 1 

1 J = 1 1 

DO 5 J= I , N1 

Z\ IJ) = DOT2(NN2,XI I ) »Y(J) ) 

5 I J= I J *-ND I M 
J= I I 
I J= J 

3 U=IJ-NDIM 

IFIIJ.LT. I) GO TO 10 

J = J-1 

zcr j ) = z ( j) 

GO TO 3 
10 11 = 1 I +ND IM 1 
RETURN 
END 

FUNCTION QOT2(NN,A,B) 

DIMENSION A{ 1 } ,91 1 j 
COMMON/MA I N 1 /NDI M 
QOT2=0 

DO l 1= l , NN , NO I M 
OOT2=nOT2+A( I ) *9 ( I ) 

1 CONTINUE 
R FTURN 
END 

SUBROUTINE REGSIM 

SUBROUTINE REGS I M< N, M, IR , ACL , G, C , XO , DT ,NS ,NP ) 

DIMENSION ACL(l) ,G( 1 > ,C(1) ,X0( 1) 

C0MMDN/MAINI/NDIM,DUM1 ( 1 ) 

COM M ON/ I NOU/K IN,KOUT 
IFINDIM.LT. N) WRITEIKOUT, 1000) 

IF(NDIM.LT.M) WRI TE I KOUT , 1000 ) 

IFINDIM.LT. IR) WRITEIKOUT, 1000) 

IFINDIM.LT. M) CALL EXIT 
IFINDIM.LT. M) CALL EXIT 
IFINDIM.LT. IR) CALL EXIT 
1000 FORMAT!/, 16H DIMENSION ERROR,/) 

LOGICAL CONT 
C JNT= . TRUE. 

CALL LNS I M2 IN, IR , ACL , C , XO, DT ,NS ,NP , M ,G ,CONT ) 

RETURN 

END 

SUBROUTINE LNS I M2 

SUBROUTINE LNS I M2 ( N, R , A, C , X 0, DT , NSTEPS , NPRPL , M, G , CONT ) 
DIMENSION A(l) ,C( 1) ,X0( 1) ,G(1) 
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DIMENSION X(30) » Y ( 30) , INDEX (30 ) ,NSYM( l ) ,U(30) 
INTEGER R »RPl 
LOGICAL CONT 
COMMON/INOU/KINfKOUT 
COMMON/PLOT/SCALE»ARRAY( I) 
C0MM0N/MAINI/NDIM,DUM1( I) 

C0MM0N/MAIN2/EAU) 

LOGICAL PRINT, PLOT 
I F (CONT ) CALL MSCALE ( M, N » G ,-l . 0 ,G ) 

RP1=R+1 
NR*NSTEPS+i 
NC=R+ 1 

I F ( CONT } NC=NC+M 
PLOT= .TRUE. 

PRINT*. TRUE. 

IF(NPRPL.LT.O) PRINT*. FALSE. 

IF(NPRPl.FQ.O) PLOT*. FALSE. 

IF (NPRPL.EQ.O) GO TO I 

I F ( N STFPS . GT. { 50*SCALE)) WRI TE ( KOUT ,2000) 

IF( N STEPS. GT.t 50*SCALE)) CALL EXIT 
2000 FORMAT!/, 12H SCALE ERROR,/) 

I CALL MEXP (N , A ,0T, EA ) 

T = 0. 

CALL VMMUL ( R , N,C , XO , Y ) 
tF(CONT) CALL VMMUL ( M, N, G , XO , U) 

I F ( .NOT. PRINT) GO TO 5 
I F ( .NOT. CONT) WRITE (KOUT, 80 ) 

I F ( CONT ) WR I TE ( KOUT , 85 ) 

W R I T F ( KOUT , 90 ) 

I F ( CONT ) WR I TE ( KOUT ,95 ) 

WRITE(KOUT, 100) T 

WRITE (KOUT, 105) ( Y( I ) ? 1= L »R ) 

T F (CONT ) WRITE (KOUT, 103) ( UC I), 1=1, Ml 

5 I F ( .NOT. PLOT) GO TO 9 
AP. RAY(1)=T 

DO 6 J=2,RP1 

6 ARRAY ( l-M J— 1 ) *NR) = Y( J-l ) 

I F ( -NOT. CONT) GO TO 9 
DO 7 J"- 1 , M 

7 ARRAY ( 1+{R+J)*NR)=U(J) 

9 00 60 K=1,NST5PS 

CALL VMMUL(N,N,EA,X0,X) 

CALL VECTEQ! N,X,XO) 

CALL VMMUL (R,N,C,XO,Y) 

T=K*DT 

I F ( CONT ) CALL VMMUL ( M , N , G , X , U ) 

I F ( .NOT. PRINT ) GO TO 54 
52 WRITE(KOUT, 100) T 

WRITEfKOUT, 105) ( Y ( I ) ,I=1,R) 

IF (CONT) WRITE (KOUT , 105 ) ( U( I) , I = 1 , M) 

54 l F { .NOT. PLOT) GO TO 60 
DO S3 J*2,RPl 

^8 ARRAYl 1 + KM J-l )*NR)=Y( J- U 
ARRAY ( 1+K)=K*CT 
I F ( .NOT. CONT) GO TO 60 
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DO 59 j* i , M 

59 ARRAY! 1+K+! R+J ) *NR ) =U ( J) 

60 CONTINUE 

IF{C0NT) CALL MSCALEtM, M, 0,-1. 0.6) 

I F ( - NOT. PLOT) GO TO 120 

NSYM! l )= 13 

NVAR=1 

NH0RZ=50*SCALE-H 
NGRI0H=5*SCALE 
NSOR T=0 
DO 65 1 = 1, R 

indexii)=i 

WRITE!KOUT, 110) I 

65 CALL PRNPLT ( NR, NC.NVAR, ARRAY? INDEX, NHORZ,NGRI DH, NSORT, NS YM) 
I F ( . NOT . CONT) GO TO 120 
NSYM! 1 ) = l 2 

00 70 1=1, M 

1 NOE X ( l )=R*I 
bRITEtKOUT, 115) I 

70 CALL PRNPLT ( NR, NC,NVAR, ARRAY, INDEX, NHORZ , NGRI OH, NSOR T , NSYM ) 
30 FORMAT!/, 28H SIMULATION OF LINEAR SYSTEM,/) 

85 FORMAT! /,31H SIMULATION OF LINEAR REGULATOR,/) 

90 F0RMATU1H OUTPUT Y) 

95 FORMAT! 1ZH CONTROL U) 

100 FORMAT!/, 5H T = ,F5.?,/) 

105 FORMAT! I X ,0P1 0E13.4) 

110 FORMAT (1H1,/ ,31X»IHY» I2,12H VERSUS TIME) 

115 FORMAT! IH1,/,31X, IHU, 12, 12H VERSUS TIME) 

120 CONTINUF 

RETURN 
END 
C 

C SUBROUTINE MEXP 

C 

SU3R0UTI NE MEXP ! N, A, T, EA ) 

DIMENSION A(i) ,EA( 1 ) , C ! 30 ) , D ( 3 1 ) , E I 30) 

COM MQN/MAI Nl/NDIMfX(l ) 

NO IM 1=N0 I M+ 1 
NN=NO I M*N 
N M 1 = N- 1 

IF (N.GT . 1) GO TO 5 
EA! I )=EXP( T*A( 1) ) 

RETURN 
5 W= 1 • 0 

DO 10 1=1 ,NN , NDI M 
I L = I +NMI 

00 10 J=I , IL 
10 EA! J )-A\ J ) 

C 1=XN0RM( N , A ) 

1 ND=0 
L= 1 

T 1 = T 

15 IF(4BS!T1*C1) .LE.3.0) GO TO 20 
T l=T 1/2 . 

I N0= I MD+ 1 
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GO TO 15 
20 C2*0. 

DO 25 !-L v NNfNDIML 
25 C2=C2-EA( I ) 

C 2*C 2/FLOAT ( L ) 

C ( L ) *C2 
0(L*1)=0. 

1 I =N + L-L 

e( m*w 

I 1 3 L 

00 35 I*l#NNfNOIM 
I L=I+NMl 

DO 30 J=I f 1L 
30 X(J)=FA(J) 

Xt m = X( II)+C2 
35 l 1 = 1 H-NDIM1 

IF(L.EQ.N) GO TO AO 
CALL MMUL( X,A»N,N,N»EA) 

M*W*T l/FLQAT ( L ) 

L = L + l 
GO TO 20 
AO CONTINUE 

C CAN CHECK X:0 FOR ACCURACY 

J = N+ 25 
DO 50 L = N,J 
DO A5 K= 1 * N 

n C K >* ( Of K*>1 )-W*C(K) )*T L/FLOATtL ) 
A5 E(K)=E(K)+D(K) 

50 W=0( 1) 

1 1 = 1 

00 AO I=l,NN,NDIM 

1 L = I +-NM 1 

DO 55 J=l » IL 
55 EAIJ)=E( 1)*A( J) 

E A ( 1 1 ) = EA ( I I ) + E ( 2 ) 

60 I I = I I *ND I Ml 

IF(N.FQ.Z) GO TO 85 
DO SO L= 3 » N 

CALL MMULIEA, AtNtN,NfX) 

I 1 = 1 
C 2 = E ( L ) 

00 75 I«1,NN»NDIM 

1 L = I +NM 1 

DO 70 J=I , IL 
70 5 A ( J ) = X ( J ) 

E A t I I ) = EA ( I I ) +C2 
75 1 1 =1 1 +NOI M 1 
80 CONTINUE 

S5 I F { IN0.5Q.0) RETURN 

00 100 L= 1 f INO 

DO 90 I = l t NN, ND I H 

1 L = I *-NM 1 

00 90 J=I , IL 
90 X{ J ) =FA ( J ) 

100 CALL MMUL(X,X,N,N,N,EA) 
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RETURN 

END 

SUBROUTINE PRNPLT 

SUBROUTINE PRNPLT ( NR , NC, N VAR » ARRAY, I NDEX, NHORZ » NGRI DH ,NSORT » /VS 
DIMENSION ARRAY! L) , I NOEX { l) , KAXI SC 101 ) ,NSYM( I ) , YLABEL ( 1 1 ) 
INTEGER CHAR! 15) 

COMMON/ I NOU/K IN* KOUT 

DATA CHAR! I ) « CHAR ( 2 ) » CHAR ( 3 ) » CHAR ! A) , CHAR! 5) , 

*CHAR(6),CHAR!7) , CHAR! 8) , CHAR! 9) , CHAR (10) ,CHAR{ LI) , 

♦ CHAR! 12) t CHAR! 13), CHAR! 14), CHARI 15 ) / 1HI , IH2 , IH3, 1H4, 

*1H5» IH6, 1H7, IH8 , IH9, LHX, 1HE , 1HU, 1HY, 1HR , LH / 

DATA ISTAR, IPER , l DASH, IBLANK/IH* , 1H. , IH- , 1H / 

XMIN=ARRAY( L) 

XM AX=ARRAY ( NR ) 

IF(NSORT.EQ.O) GO TO 10 
DO 5 1=1, NR 

I F! XMIN.GT .ARRAY ! 1 ) ) XMI N=ARRAY ! I ) 

5 IFfXMAX.LT. ARRAY! I ) ) X MAX=ARRAY ( I ) 

10 Y^AX= ARRAY ( l + INDEX ( 1 ) *NR ) 

YMIN=YMAX 
DO 50 J= 1 , NVAR 
DO 50 1=1, NR 
I J=IMNDEX( J)*NR 

IFIYMIN.GT .ARRAY! I J ) ) YM I N=ARRA Y! I J ) 

50 IF! YMAX.LT. ARRAY! I J) ) YMAX=ARRAY ( l J ) 

IF(A3S!YMAX-YMIN) .GE.1.0E-20) GO TO 55 
YMAX=YMAX+ l . 

YMIN=YMIN-L. 

55 CALL RNDOFF ( YMAX, YMIN ) 

YLABEL ( 1 )=YMIN 
YLABf- L! 6)=YMAX 
DO fcJ 1=2,5 

60 YLABEL (I )=( 1-1 )*{YMAX-YMIN) /5.0+YMIN 
WR ITEUOUT, 500) (YLABEL! t ), 1=1,6) 

I F ! NGRI DH. G6 r L) GO TO 65 
L INEPR=NH0RZ/5 
GO TO 70 

65 L INEP R=NGRI DH 
70 KL INE=1 
L I NE= 1 

S TEP X=l XMAX-XMIN) /( NHORZ— 1 ) 

STEP Y= ( YMAX— YMI N ) / 50. 

DO 220 I ND= 1 , NR 

K STEP=! ARRAY ( I NO ) -XMI N) /STEPX+l . 5 
75 IF (LINE-KL1NE) 140,90,90 
30 XLA9FL=STEPX*!KLINE— 1) +-XMIN 

l F! (NGRIDH.EQ.O) .AND. (LINE.NE.l) ) GO TO 130 
DO 90 1=2,51 
90 KAX IS! I )=I DASH 
DO 100 1=1,51,10 
100 < AXIS! I ) = riTAR 

105 IF(KSTEP-LINE) 110,110,125 
1 10 on 120 J= i « NV AR 
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K* ( ARRAY ( I NO* INDEX { J )*NR ) -YM IN ) /STEPY+l . 5 
120 KAX I S (K ) *CHAR ( NSYM ( J ) ) 

125 WR ITF ( KQIJT » 600 ) XLABEL, { KAX I S( I ) , 1*1 , 51 ) 
KLINE*KUNEfUNEPR 
GO TO 200 
130 OH 135 1*2, 51 
135 KAXIStI)* I BLANK 
K A X I S ( 1 )*IP5R 
GO TO 105 
1 AO 00 150 1*2, 5L 
150 KAXISU l«IRLANK 
KAXIS( 1)*IRER 
IP(NGRIDH.FQ.O) GO TO 165 

00 160 1*11,51,10 
160 KAXI S( I ) = IPER 

165 IF(KSTEP-LINE) 170,170,190 
170 00 130 J*i,NVAR 

K= { ARRAY ( INDUNDEXt J)*NR)-YMIN)/STEPYU.5 
180 K AXI S (K )*CHAR ( NSYM ( J ) ) 

190 rtR IT E ( KOUT ,700) (KAXISU ), 1*1,51 ) 

200 L I N£=LI NF + 1 

1 F (L JNE-NHORZ-l ) 210,210,230 
210 IF(KSTEP-LINF) 22 0,75,75 

220 CONTINUE 
230 RETURN 

500 FCT M AT(//,9X,6( 1PE 1 0* 2 ) ) 

600 r ORMAT ( IPE13.2,2X»51A1) 

700 FORMAT ( 15X, 51 A1 ) 

RETURN 

END 

SURRQUT INF MSCALE 

SUBROUTINE MSGAL E ( Ml , N2 , A , X , B ) 

DIMENSION A( 1) ,B( 1) 

CO M MON/ M A I Nl/ND IM 
JEND=N2*NDIM 
DO 10 1*1, N1 
00 10 I J=I , JEND,NOIM 
10 B ( I J ) =X*A ( I J ) 

RETURN 

FNO 

SUBROUTINE VMMUL 

SUBROUTINE VMMUL (Nl ,N2 ,X,Y, Z) 

DIMENSION X(l),Y(l),Z(l) 

COMMON/M A IN 1 /NDI M 
DO 10 1=1, Nl 
Z ( I ) =9 

DO 10 j - 1 ,N2 

10 Zl l )=Z( I )+X( l+( J-l)*NDlM)*Y(J) 

RETURN 

END 


95 


o o n o o o 


Appendix B LINEAR Computer Program 


C SUBROUTINE VECTEQ 

C 

SUBROUTINE VECTEOI N, X, Y) 

DIMENSION XU) ,Y( L) 

DO LO I* L , N 

10 Ym*xm 

RETURN 
END 

SUBROUTINE RNDOFF 

SUBROUTINE RNDOFF ( DMAX, DM IN ) 

COMMON/ INOU/K I N.KOUT 
RANGE*DMAX-DMIN 
J=-l 
K=0 

IFIRANGE.LE.2. ) J = L 
DO LO 1=1,25 

IF{ (RANGE. GT. 2. ) .AND. (RANGE. LE. 20. ) ) GO TO 20 
K=K + J 

10 R ANGE=RANGE*IQ.**J 
W Rl T F ( KOUT , 100) 

CALL EXIT 
20 T=DMIN*LO.**K 

I F(DMIN.LT.O) GO TO 30 
N=T 

DMIN=N*10.**(-K) 

GO TO 40 
30 N = T 
T L = N 

IF(Tl.GT.T) N=N-l 
DMIN=N*10.**(-K) 

40 T=QMAX*10.**K 

I F ( DMAX . GT . 0 ) GO TO 50 
N=T 

,DMAX=N*10.**(-K) 

GO TO 60 
50 N=T 
Tl = N 

IF(Tl.LT.T) N=N+L 
OMAX=NnO.**(-K) 

60 CONTINUE 

100 FORMAT (/fix, 32HNUMBERS CUT OF RANGE FOR PLOTTER,/) 
RETURN 
END 

SUBROUTINE CBS 

SUBROUTINE QBS ( N, l R , A, C , NOS ) 

DIMENSION A ( 1 ) , C ( 1 ) 

COMMON/ I NOU/K IN , KOUT 
C0MM0N/M4IN1/NQIM 
I F(NDIM.LT.N) WRITE (KOUT, 1000) 

IF(NOIM.LT. IR) WRITE (KOUT, 1000) 

IF( NDIM.LT.N) CALL FXIT 
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I FINDIM.LT. IR) CALL EXIT 
1000 FORMAT ( / , 16H DIMENSION ERRORf/) 

CALL TRANSK N» A, A) 

CALL TRANS3I IR,N,C,C) 

CALL CONT(N,lR,A,C,NOS) 

IFINOS.LT. N» GO TO 10 
WRITEIKOUT, 20) 

GO TO 50 

10 WRITEIKOUT, 30) 

WRITEIKOUT, 40) NOS 

20 F0RMAT(/,10X, 20HSYSTEM IS OBSERVABLE , / ) 

30 FORMAT ( / , LOX ,22HSYSTEM IS UNOBSERVABLE,/) 

40 FORMAT ( 15X, 30HNUMBER OF OBSERVABLE STATES = ,12,/) 

*0 CALL TRANS1(N,A» A) 

CALL TRANS3(N,IR,C»C) 

RETURN 
END 

SUBROUTINE TRANS 1 
SUBROUTINF TRANSl ( N* A , AT ) 

SETS AT= ATRANSPOSE A AND AT CAN BE EQUIVALENT, BOTH ARE SQUARE 
DIMENSION A ( 1 ) , ATI I ) 

C3MM0N/MAINI/NDIM 
ND I.M 1=ND I M+l 
NN=N*NDIM 
DO I 1=1 ,NN ,NDI Ml 
I J=I 

DO 1 J 1= I ,NN, NDIM 
TEMP=AI I J) 

AT ( I J ) = A ( J I ) 

ATI JI )=TEMP 
IJ=I J+l 
l CONTINUE 
RETURN 
END 

SUBROUTINE TRAN S3 

SUBROUTINE TR ANS3I NR, NC , A , A T ) 

C SETS AT= TRANSPOSE OF A , A AND AT CAN BE EQUIVALENT , A IS NR * /VC 
DIMENSION All ) ,AT(1 J 
COMMON/MAIN 1/ND I M 
ND IM1=NDI 5 + 1 
NR L=NR- 1 
NC1=NC-1 

IFINR.GE.NO N=NC 
IFINR.LT. NC) N=NR 
NN=N*NDIM 
DO 1 I=1,NN,NDIM1 
I J= I 

DO 1 JI=I,NN» NO I M 
T CM P'=A ( I J) 

A T 1 1 J ) = A ( J I ) 

AT(JI)=TFMP 


97 


Appendix B LINEAR Computer Program 


I J^I J+l 

1 CONTINUE 
IF(NR-NC) 5,2,3 

2 RETURN 

3 DO 4 I=NC , NR1 
DO 4 J= 1 »NC 

4 A( J + I*ND IM 1 = A C I H*( J-l) *NOIM) 

RETURN 

5 DO 6 1*1, NR 
DO 6 J=NR , NC I 

6 AC J+l+< I-l )*NDIM)=A( H-J*NDIM) 

RETURN 

END 

C 

C SUBROUTINE CONT 

C 

SUBROUTINE CONT ( N» M, A, B *NCS ) 

DIMENSION AC l ) , B( 1 ) 

COMMON/ MA IN L/NDIM, DUML (1) 

CALL EQUATE (N,M»DUMl » B ) 

20 INOEX^L 
I R - M 

CALL ORTHNM ( DUM I , N , IR , 0 ) 

NC 5= I R 

IF (IR.EQ.N) RETURN 
I RL AST = 0 
30 INDEX = INDEX+L 
ML = I RLAST+l 
I RL AST = IR 

C MULTIPLY COLUMNS 'ML' THROUGH ' IRLAST ' 3Y MATRIX A 

K L=< Ml-L )*NDI M+L 
K2=IR*NDIM+L 

3 L ML = MINOC IRL AST+L-M 1 , N- IR > 

33 CALL MMUL I A, DUMl ( KL ) , N , N ,ML , OUM l ( K2) ) 

35 ML = ML+ML 
I WOLD = IR 
IR = IR+ML 

CALL ORTHNM ( DUM L , N, IR , I R HOLD ) 

NC S= I R 

IF (IR.EQ.N) RETURN 
IF (ML-IRLAST) 31,31,50 
50 IF (IR-IRLAST) 60,60,30 
6 INDEX = 0 
RETURN 
END 
C 

C SUBROUTINE ORTHNM 

C 

SU3R0UT INE ORTHNM ( A, NR, NC ,NNCO ) 

DIMENSION ACL) 

C3MM0N/MAINL/NDIM 

C COLUMN ORTHO-NORMALIZATION ROUTINE { GRAM- SC HM I DT ALGORITHM) 

C NR IS THE COLUMN LENGTH 

C NC IS THE NUMBER OF COLUMNS 

C NNC 0 IS THE NUM3ER OF ORTHONORMAL COLUMNS AT CALL TIME 
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C AND: NC WILL BE REDUCED IF LINEARLY DEPENDANT COLUMNS OCCUR. 

NNC=NNC0 
l II 3 1 

A IF ( NNC ) 5,21,5 
5 IF ( NC-NNC ) 50,50,10 
10 ILIM=NDIM*NNC 
II = I L I M 1 
15 DO 20 K* 1 , 2 

DO 20 I=1,ILIM,NDIM 
W=-D0T(NR,A( I ) ,A( II) ) 

CALL WAOO ( A ( II), A( I ) , A ( II) ,W,NR) 

20 CONTINUE 

21 W=D0T{^R, A( II ) ,A( I I ) ) 

IF (W - 1. 06-12) 30,40,40 

30 IF (NO-NNC-L) 33,33,35 
33 NC = NNC 
RETURN 

35 I — t NC— 1 ) *ND I M+l 

CALL VECTEQINR, A( I ) ,A( II) > 

NC-NC-1 
GO TO 4 

40 W = 1 . 0/ SORT ( W ) 

CALL MSCALEINR, i,A(II) ,W,A( I I) ) 

NNC=NNC+ l 
GO TO 5 
50 RETURN 
FND 

SUBROUTINE WADD 

SUBROUTINE W ADD< A , B , C, W, N ) 

DIMENSION A(l),B(l),C(l) 

DO 10 1=1, N 
10 C ( I )=A ( I ) +W*B ( I ) 

RETURN 
END 

SUBROUTINE VPCTIO 

SUBROUTINE VECT IOC N, X, 10 ) 

DIMENSION XI l) 

CD M MON/INOU/KIN , KOUT 
IFII0.EQ.4) GO TO 40 
l F ( I 0. EQ . 3 ) GO TO 20 
5 CONTINUE 

INPUT 

RCAO (KIN, 1000) (XII) ,I=i,N) 

I F C 10. EQ. I ) RETURN 

OUTPUT 

20 « R I T5 ( ROUT , 1001 ) I X ( I ) , 1 = l ,N ) 

GO TO 50 
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I 


TITLE 

40 READtKIN, 1002) 

WRITE! KOUT • 1003) 

WRITEtKOUT, 1002) 

WRITEtKOUT, 1004) 

GO TO 5 
50 CONTINUE 

1000 FORMAT l BE 10 • 0 ) 

1001 F0RMAT(IX,1P10E13.4) 

100? FORMAT! IX, 79H 

C ) 

1003 FORMAT!//) 

1004 FORMAT!/) 

RETURN 

END 

SUBROUTINE CON 

SUBROUTINE CON ( N, M, A, B , NC S) 

01 MENS ION A! I ) ,B( 1) 

COMMON/INOU/KIN.KOUT 
COM.MON/MA IN1/N0IM 
I FINOIM.LT.N) WRITEtKOUT, 1000) 

IF(NOIM.LT.M) WRITEtKOUT, 1000) 

IF(NOIM.LT.N) CALL EXIT 
I F( NOIM.LT.M) CALL EXIT 
1000 FORMAT ( / , 16H DIMENSION ERROR,/) 

CALL C0NTIN»M»A»8»NCS) 

IFINCS.LT.N) GO TO 10 
WRITEtKOUT , 20 ) 

RETURN 

10 WRITE! KOUT , 30 ) 

WRITEtKOUT, 40) NCS 

20 FORMAT!/, 10X , 22HSYSTEM IS CONTROLLABLE,/) 

30 FORMAT!/, L0X,24HSYST6M IS UNCONTROLLABLE ,/ ) 

40 FORMAT! I 5X , 32HMUMBER OF CONTROLLABLE STATES = ,12,/) 
RETURN 
END 

SUBROUTINE TRANS2 

SUBROUTINE TRANS2 ( NR,NC, A , AT > 

C AT=ATRANSPOSE, A IS NR X NC 

DIMENSION A! 1) , AT ! 1) 

COMMON/MAINl/NDIM 
NN-NC*NO I M 
I 1 = 1 

DO 6 1=1, NR 
J J = I I 

DO 5 J=I , N N , N D I M 
AT! JJ) = A( J) 

5 JJ=JJ+l 
I I=t I+NDIM 


i 
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6 CONTINUE 
RETURN 
END 

SUBROUTINE EIGVAL 

SUBROUTINE EIGVAL(N,A) 

DIMENSION AU),RR(30),RI(30) 

COMMON/MAIN I /ND I M 
C OMMQN/ I NOU/ KIN, KOUT 

I F(NDIM.LT.N) WRITEIKOUTtlCOO) 

IF(NDIM.LT.N) CALL EXIT 

1000 FORMAT { / » 16H DIMENSION ERROR,/) 

CALL MEIGVIN, A,RR,RI ) 

WRITEIKOUT, 100) 

WRITE (KOUT, 200) (RR(I),RI(I),I-l,N) 

100 FORMAT ( / , I9H MATRIX EIGENVALUES, / ) 

200 FORMAT ( I3H REAL PART = ,E10.3,13H l MAG PART = ,E10.3) 
RETURN 
END 

SUBROUTINE VMAT2 

SUBROUTINE VM AT2 ( Nl , N2 , N3, A , B , C , 0, E ) 

DIME NS I ON A(L),B(1),C(L),D(1),E(L) 

COMMON /MAIN l/NDIM 
00 30 1=1, NI 
TEMP=0. 

DO 10 J=1,N2 

10 TEM=> = TEMP+A( I+I J-l)*NOlM)*S( J) 

DO 20 J=i,N3 

20 TEMP=TEMP+C( I+( J-l )*NDIM)*D( J) 

30 F ( I ) =TEMP 
RFTIJRN 
END 

SUBROUTINE MATS 

SUBROUTINE MAT 5 (Nl,N2»N3»X»Y,Z) 

DIMENSION XI 1) ,Y(I) ,Z( l) 

COMPUTE Z=X'Y WHERE X IS NLXN2, Y IS NIXN3,Z IS N2XN3 
CO M MON/MA I N 1 /ND I M 
DO 2 1=1, N2 
DO 2 J=1,N3 
TM = 0. 

II = ( 1-1 )*NDIM 
JJ = ( J-l)*NDIM 
^0 1 K= 1 , N 1 

1 TM = T M+X ( I I+K ) *Y( K + J J ) 

2 Z ( I J J ) = TM 
RFTURN 
END 

SUBROUTINE MAT3A 
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SUBROUTINE MAT3AIN1,N2, X, Y, Z ) 

C Z=X*YX Y=Y * IS N2XN2 X* IS NIXN2 

DIMENSION X(I),Y( 1) ,Z( l) 

COMMON/MAI Nl/NDIM 
NDIM1=NDIM+1 
NN1=N1*NDIM 

CALL MMUL(Y,X,N2,N2,Nl, Z J 
1=0 

DO LO 1 1 = I , NNi » ND IM 
J=II+l 
I J=J-I 

DO 5 JJ=I 1 1 NNI * NOIM 
Z( J)=D0TIN2,XI 1 1 ) ,Z(JJ) ) 

5 J=J+I 
1=1 + 1 
J J = I 

DO 10 J=II , IJ 
Z( J)=Z(JJ) 

10 J J=J J+ND I M 
RETURN 
END 

SUBROUTINE LINSIM 

SUBROUTINE L INS I M ( N , I R, A, C , XO, DT, NS ,NP ) 
DIMENSION ACL) ,CU) ,X0ll) ,Gll) 

COMMON/MAIN l/ND I M , DUMl ( 1 ) 

COMMON/ INOU/KINi 'TOUT 
IFINOIM.LT. N) WRITEIKOUT, 1000) 

IFINDIM.LT. IP.) WRITEIKOUT, LOOO) 

IFINOIM.LT. N) CALL EXIT 
IFINDIM.LT. IR) CALL EXIT 
1000 FORMAT!/, 16H DIMENSION ERROR,/) 

LOGICAL CONT 
CONT=. FALSE. 

M=l 

CALL LNSTM2(N,IR,A,C,XO,DT,NS,NP,M,G,CONT) 
RETURN 
END 

SUBROUTINE EQCOST 

SUBROUTINE EQCO ST ( N , M, A, 3 , R , S,Q,AEQ,QEQ) 
DIMENSION M l ) ,Btl ) ,RI l) , St 1) ,Q( 1) 

01 MENSION AEQ I l ) ,QEQ( L ) 

COMMON / I NOU/.K I N , KOUT 
COMMON/MAINl/NDI M, DUM1 1 1 ) 
COMMON/MAIN2/DUM2 I 1 1 
1 FINDIM.LT.N) WRITEIKOUT, 1000) 

IFINOIM.LT. M) WRITEIKOUT, 1000) 

IFINDIM.LT. N) CALL EXIT 
IFINOIM.LT. M) CALL EXIT 
1000 FORMAT!/, I6H DIMENSION frrqr,/) 

CALL EQUATE l M »M»DUM2 » R ) 

CALL GM I NV I M , M , DU M2 , OUM 1 , MR , 0 ) 
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IFIMR.EQ.M) GO TO 10 

wRiretKouT, ioo) mr 

CALL EXIT 

10 CALL MAT4 ( N , M» DUMl , S» 0UM2 ) 

DO 20 I* 1 , N 
00 20 J* 1 »N 

20 3 BO ! I* ( J- l ) *ND I *1 )«Q ( 1 + ( J- L ) *ND I H ) -DUM2 I I -M J- 1 ) ♦NO I M ) 

CALL HAT6(M,M,N,DUML f S.0UM2) 

CALL MMUL ( 3 , DUM2 , N» M, N» DUMl ) 

DO 30 1*1, N 
00 30 J*l,N 

30 ASQI I-M J 1 )*NDIM!*A(lMJ-l)*NDIM)-DUMlUMJ-l)*NDIM) 

100 FORMAT (/ » 25H ERROR MATRIX R HAS RANK ,I2,12H LESS THAN M) 
RETURN 
FNO 
C 
C 
C 

SUBROUTINE TITLE 
COMMON/ I NOU/K IN, KOUT 
READ (KIN, LQ02 ) 

W R I T E ( KOUT , L003 ) 

WRITEIKOUT, 1002) 

WRITE (KOUT, 1004) 

1002 FORMAT! IX, 79H 

C ) 

1003 FORMAT ( // ) 

1004 FORMAT!/) 

RETURN 

END 
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FILE INTODEl FORTRAN AL 

PROGRAM TO INTEGRATE GENERAL ODE'S 

THIS VERSION WAS USED ESPECIALLY TO GENERATE TRAJECTORIES FOR 
THF LANGLEY/NFWSOM CONTRACT SUMMER 1980 

D IMF NS ION XXI 2010) »YYI 2010) vNPTSI 10) ,XYl (4) 

DIMFNSION NXPLI 10) ,NXPR( 101 , X0( 12 ) , X I ( 12 ) 

DIMENSION X(l2)fF!10)t XNA ME S ( 10 ) » XL 9L S ( 5 » 10 ) 

DIMENSION CNAMESI 10) ,XPLBL( 10) , YPLBL ( 10 ) , T1 TL I 20 ) ,NCINDX{ 10) 

DATA ISTART/1/ 

DATA XYL/O. ,5.0,-l.,l./ 

OAT A WIDTH,FEIGHT,TICKL/7., 5-,. 08/ 

OAT A NXPR/l,2,3,4,5,5*0/,NXPL/l,2,3,7*0/ 

DATA NSTEPStNPLOTf NPRI NT/3* 1/ 

COMMON/CASE/NCASE 

COMMON/NNOUT/NOneV 

C0MM0N/MAT/A(4,4) ,B(4,3),XK( 3,4),NX,NU,NINT,NT0T 

NAMFLIST/INTODE/XI.TMAX, I PR I NT , TP<U NT ,DT , A, B , XK, NX, NU , N INT , N TOT, 
1 NXPR,IP10T.TPL0T,N0DEV, I CONT,NXPL, WIDTH, HEIGHT, TICKL, 
l NCHARS , NL INES , X NAMES ,XL BLS , XPL BL » NXC, YPLBL , NYC , T I TL ,NTTL 

INPUT VARIABLES IN NAMELIST: 

XII ) = INITIAL VALUES FOR STATE VARIABLES 

TMAX = MAXIMUM (FINAL) INTEGRATION TIME 

IPRINT = INDEX TO DETERMINE THE AMOUNT OF PRINTOUT 

1, NORMAL OUTPUT 
0, MINIMUM OUTPUT 

2, EXTRA OUTPUT 

TPRINT = TIME INCREMENT FOR TRAJECTORY PRINTOUT IN TABLE 
DT = INTEGRATION STEP SIZE, SHOULD BE A WHOLE DIVISOR 
OF TPRINT AND TPLQT 
A( , ) = NX BY NX SYSTEM MATRIX 
B( , ) = NX BY NU CONTROL MATRIX 
XK{ , ) = NU BY NX FEEDBACK MATRIX 
NX = NO. OF STATE VARIABLES 
NU = NO. OF CONTROL VARIABLES 

NINT = NUMBER OF VARIABLES TO BE INTEGRATED IN RKINT 
NTOT = TOTAL NUMBER OF STATE, CONTROL AND AUXILLARY VARIABLES 
NXPRI ) = ARRAY OF INDICFS OF VARIABLES IN TRAJECTORY TABLE 
I PLOT = INDEX THAT CONTROLS THE GENERATION OF PLOTS 

0, NO PLOTTER OUTPUT 

L, GENERATE PRINTER PLOT ONLY 

2, GENERATE BOTH PRINTER PLOT AND VERSATEC PLOT 
TPLOT = TIME INCREMENT FOR STORING POINTS FOR PLOTTER 
NOPFV = UNIT NUMBER OF OUTPUT DEVICE, NOMINALLY 6 
ICONT = INDEX TO CONTINUE OR STOP MULTIPLE CASE RUNS 

3, STOP AFTER THIS CASE 

1, CONTINUE TO NEXT CASE 

NXPLI ) = ARRAY CF INDICES OF VARIABLES TO BF PLOTTED 
WIDTH = WIDTH OF VERSATEC PLOT IN INCHES 
HEIGHT = HEIGHT OF VERSATEC PLOT IN INCHES 

TICKl = LENGTH OF TICK MARKS AND HEIGHT OF CHARACTERS (INCHES) 

IF NEGATIVE, PLACE TICK MARKS ON INWARD SIDE OF AXIS 
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C IF POSITIVE. PLACE TICK MARKS ON OUTSIDE OF AXIS 

C NCHARS * NO. OF CHARACTERS PER LINE IN PRINTOUT (WIDTH) 

C NLIN6S * NO. OF LINES FOR PRINTER PLOT (HEIGHT) 

C XNAMFS ( ) * ARRAY OF ONE-WORD (FOUR-LETTER) NAMES OF VARIABLES 
C XL9LS( 5 » ) * ARRAY OF 5-WORD (20-CHARACTER) DESCRIPTIONS OF VARIABLES 
C XPLBL ( ) * ARRAY CONTAINING THE X-AXIS LABEL FOR PLOTS 
C NXC * NO. OF CHARACTERS (LETTERS ) IN X-AXIS LABEL 
C YPLBL ( ) * ARRAY CONTAINING THE Y-AXIS LABEL FOR PLOTS 
C NYC * NO. OF CHARACTERS IN Y-AXIS LABEL 

C T ITL ( ) * ARRAY CONTAINING A TITLE TO APPEAR IN PRINTOUT ANO ON PLOTS 

C NTTL * NO. OF CHARACTERS IN TITLE 
C 

NPMAX * 2010 
NODEV * 6 
NCASF * 0 

1 READ(7, INTOOE) 

IF( I PR INT .GE.2)WRITE(N00EV» INTODE) 

NCASE * NCASE *• 1 

DETERMINE NO. OF PRINT COLJMNS 

DO 2 1=1,10 

I F ( NXPR ( I ) .LE.OJGO TO 3 
CONTINUE 
NPRTOT = 1=1 

DFTPRMINF NO. OF CURVES ON PLOT 

DO 4 I = 1,10 
I F( NX°L ( I ) .LE.OJGO TO 5 
CONTINUE 
NCURV = I-i 

DFTPRMINF NO. OF INTEGRATION/PLOT/PRINT STEPS 

IF(DT.NE.O. )G0 TO 6 
WRITE(NODEV» 104) 

FORMAT ( / , • *** FATAL ERROR *** DT IS ZERO. GO TO NEXT CAS£.‘) 
GO TO 30 

IF( IPLOT.LE.OJTPLOT * TPRINT 
TPLOT = AMAX 1 ( TPLOT , DT ) 

TPRINT = AMAXK TPRINT, TPLOT) 

NSTFPS = (TPLOT * .OOOOD/DT 
NPLOT = (TPRINT + .0000 1 ) /TPLOT 
NPKINT = ( TMAX .00001) /TPRINT 
NPNTS = NPLOT*NPRINT + 1 

SET UP DATA FOR PLOTS 

I F { NCUR V .NE . 0 ) NPNT S = MtNO( NPNTS, ( NPMAX— NCURV) / NC'JRV) 

DO 7 1=1, NCURV 
C NCINDXU) = NPNTS*( FLOAT ( I ) /FLOAT(NCURV+ l > ) 

NCINDXt I }=NPNTS/ I5.*l0*( 1-1) 

C THIS LINF PUTS THE LA8LES ON APLOT CURVES CLOSE TO THE LEFT MARW/V, 
K = NX PL ( I ) 
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CNAMESt l)»XNAMESCK) 

7 NPTS ( I ) * NPNTS 

C 

C PRINT BANNER 
C 

WRITE (NODEV, L05 INCASE 

105 FORMAT! IHl , ' ********** PROGRAM INTODE CASE* tr X 3 « / ) 

CALL DAT I ME ( I M, I D* I Y t IH* IMN»AP, ’INTODE ', NODEV) 
t F ( l PRINT. GT.O) WRITE l NODEV, LOO )TI TL , TMAX, TP R I NT » TP LOT »DT , 

1 I°R INT, IPLOT, NCURV, NPNTS 
LOO FORMAT ( / » IX ,20A4 , 

l //• TMAX , TPR INT • T PLOT , DT *•, 

l AFL2.5,//’ I PR INT » I PLOT * NCURV , NPNT5 *• ,415) 

C 

C INITIALIZE STATE* CONTROL AND AUXILLARY VARIABLES 
C 

CALL INITAL(X*XI»IPRINT) 

TIME * 0. 

CALL AUXVAR I X * T I .ME) 

IP = L 

IT! IPLOT.GT.OLCALL XYSTOR ( XX , YY , I P , T I ME, X ,NP.NT S * NPMAX , NCURV . NX PL) 
MODE s O 

CALL PT ABLE ( T IMF, X* NX PR, NPRTOT * XNA M ES» XLBLS, MOOF , NODE V, NCHARi ) 

C 

C BEGIN INTEGRATION LOOP 
C 

00 10 I* 1 , NPR INT 
DO 20 J* l , NPLQT 

00 25 K= L , NSTEPS 

25 CALL RK INTI X, TI ME ,DT ,N INT) 

CALL AUXVARt X, I IME) 

IF ( I PLOT. GT.O) CALL XYSTOR (XX ,YY , I P , T I ME , X , NPNTS , NPMAX , NCURV , /V X PL) 
20 CONTINUE 

MOOE=L 

CALL PT ABLE (T IME,X,NXPR, NPRTOT , XN AMES , XL3LS , MODE, NODE V.NCHARS) 

10 CONTINUE 

C 

C DO PLOTTING 
C 

I F { I PLOT . GE . I )CALL PPLOT ( XX , YY , NCURV* IP , XYL , NCHAR S , NL INES) 

I F ( IPLOT.LE.L ) GO TO 30 
I F ( I ST ART .GE.lJCALL PLOTS (0, 0, 50) 

IF ( I STAR T.LE.O) CALL PLOT ( WI DTH*2. , 0. ,-3) 

1 START = 0 

CALL APLOT ( XX, YY, NPTS, NCURV, XYL, WIDTH, HEIGHT, TICKL,NCASE, 
l XPLBL, NXC, YPLBL,NYC,TITL,NTTL,CNAMES,NCINDX) 

C 

C GO TO NEXT CASE IF ICONT IS GREATER THAN ZERO 
C 

30 I F ( ICONT. GE. L ) GO TO I 

IF ( I PLOT. GT.O) CALL PLOT ( 3 . ,0 . ,999 ) 

STOP 

END 

C 

C 
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SUBROUTINE RKINTIX,TIME,QT*NX) 

DIMFNSIQN XI D ,DX( 10*4) » XTEMPl 10) 

CALL XD0T(X,TIMF,0X{ l,D ) 

TIMEO * TIME 

DO 20 1*1,3 

T * DT/FLOATI 2-1/3) 

TIME * TIMEO * T 
DO 10 J* I ,NX 

10 XTRMPIJ) * X(J) ♦ T*DX I J * I ) 

20 CALL XDOT(XTEMP,TIME,DX( l,IU) ) 

DO 30 1=1, NX 

30 X(I) = X( I )*DT*(DX( I, l) + 2.*0XU , 2 )+2. *DX( I , 3) +DXI 1,4) )/6. 

TIME * TDEO ♦ DT 
RETURN 
END 
C 

SUBROUTINE XYSTOR ( XX, YY , IP , T IME , X ,NPT S, NPMAX , NCURV, NXPL ) 
DIMENSION XX U) , YY { l) ,X( l),.NXPL(l) 

00 10 I*L, NCURV 

1 I * IP •** ( 1-1 )*NPTS 

I F( I I .GT-NPMAX) RETURN 
XXIII) * TIME 
JJ * NXPL ( I ) 

10 YY(II) = X!JJ) 

IP ~ IP + i 
RETURN 
END 
C 
C 

SUBROUTINE PT ABLE ( T I ME , X , NX PR, NPR.XNAME , XL9LS, MODE* NODE V , NCHAR S) 
DIMENSION XII) * NXPR ( l) , XNAME ( I) ,X0(20) , XNAMO I 20 ) , XLBLSI B, I ) 

INPUTS: 

M00E=0, PRINT DEFINITIONS, COLUMN HEADINGS AND L ROW OF DATA 
=1, PRINT DATA ONLY 

TIME 

X{ ) = ARR AY OF STATE AND AUXILLARY VARIABLES 
NXPR ( )=ARRAY OF INDICES OF VARIABLES TO BE PRINTED 
NPR=NQ. OF VARIABLES PRINTED 

XNAMEI )=APRAY OF Y-CHARECTER LABLES OF VARIABLES 
XLBLSI 5, ) = ARRAYS OF 20-CHARACTER DEFINITIONS OF VARIABLES 
NODEV=OUTPUT UNIT NO 

NCHAR S= W I DTH OF PAPERINO OF CHARECTERS) 

COMMON/ MAT / A ( 4t 4 ) ,8(4,3) , XK I 3, 4 ) , NX , NU *NI NT , NTOT 
PRINT DEFINITIONS AND HEADINGS 
IFIMCJDE.GE.DGO TO 20 

WRITE! NODEV, 106 ) ( I, XNAME I I ) , I XL 3L S I J , I ) , J=i , 5 ) , 1=1, NTOT) 

L06 FORMAT!//,' STATF, CONTROL AND AUXILLARY VARIABLES:', 
l //,' VAR' ,4X, 'LABEL', 4X, 'DEFINITION' , 

l /, ( ' X! • ,12,' ) =' * A4, ' , • o 5 A4 ) ) 

DO 10 I = 1 , NPR 
J = NXPR! I) 
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10 XNAMOU) * XNAMEtJ) 

IFlNCMARS .IT. 130) WR I TE { NODEV , 100 ) ( XNAMOU) » 1*1 t NPR ) 

100 FORMAT (//,7X,'T IMF' ,7(6 V ,A4) , // , 1 IX, 7 ( 6X, A4) ) 

I r ! NCHARS .GE. 130 > RITEtNODEV, 105) ( XNAMO( I ) , 1*1 ,NPR) 

105 FORMAT f//,7X,»Tl^ 5 , L2(6X,A4) ,//,UX, 12(6X,A4) ) 

WR ITF(NODEV, 107) 

107 FORMAT ( 10X) 

C 

C PRINT ONE (OR MORE) ROWS OF DATA 
C 

20 NCOL* 12 

IF (NCHARS .LT. 130) NCOL-7 

NCOLl*NCOL+i 

XMIN * 10.E10 

XMAX * 10.E-10 

00 30 I -* 1 » NPR 

J * NXPrt ( I ) 

XO(U = X C J ) 

IFIXOm.NE.O.JXMIN * AMIN1 ( XMIN,ABS( X 0 ( I ) ) ) 

30 XMAX * AMAXUXMAX,ABS(XO(I) ) ) 

NPM * M I N0( NPR, NCOL ) 

C 

C OFCIDF IF OATA SHOULD 3F PRINTEO IN FLOATING POINT OR E-FORMAT 
C 

IF(XMIN. LT.. 001. OR. XMAX. GT. 99999. )G0 TO 40 
C 

C PRINT IN F10.3 FORMAT 
0 

WRITE (NOOEV,lOl)TIMC,(XO( I) ,1*1, NPM) 

101 FORMAT(IX,13F10.3) 

I F{ NPR. LE. NCOL ) RETURN 
C 

C A SECOND ROW MUST BE USED TO PRINT ALL THE DATA 
C 

wR ITE(NOOEV »l 02 )(XQ( I ) , I=NCOLL,NPR) 

102 FORMATt LIX, 12F10.3) 

RFTURN 

C 

C PRINT IN E10.3 FORMAT 
C 

40 WRITFtNODEV, 103) TIME, ( XO ( I ) , I* l ,NPM > 

LOT FORMAT! IX , F 10.3 1 1 P L 2E 10. 3) 

I F (NPP.LE.NCOL ) RETURN 
C 

C A SECONO ROW MUST «E USED TO PRINT ALL THE OATA 
C 

WR IT E ( NOQFV , 104 ) ( XQ ( I ) , I =NCOLi , NPR ) 

104 FORMAT(11X,1 3 12E10.3) 

Rr TURN 
PNO 
C 

SUBROUTINE IN ITAL ( X,X I , I PRINT) 

DIMENSION X( l ) ,Xt (1) ,DUM(4,4) 

DIMENSION 0 H 5,5) ,02( 5,5) ,03(5,6) , 04(b) ,ER( 5) ,EIl 5) 

COMPLEX EIG(5),EV£C(5,5) »C NORM 
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CQMMON/NNCUT/NODEV 

COMMON/MAT/ A( 4,4) , B( 3 ) , XK < 3, 4 ) , NX, N’J ,N I NT,NTOT 
C 

C INITIALIZE TIME AND STATE 
C 

DO ? I * ltNINT 
2 X ( 1 ) * XIII) 

TIME-O. 

CALL AUXVAR ( X » Tl ME) 

C 

C PRINT OUT A# D AND K MATRICES 
C 

CALL WRTMAT ( A ,NX,NX»NX» ' A •) 

CALL WRTMAT ( B »NX, NU, NX, • B •) 

CALL WRTMAT(XK,NU,NX,NU,*K *) 

C 

C COMPUTE CLOSED— LOOP MATRIX 
C 

DO 10 I * l , N X 
DO 10 J* l » NX 
TEMP*0 

• DO 5 K s 1 , NU 

5 TEMP-TEMP+BI I ,K)*XK(K, J) 

10 DUM(I,J)*TEMP+A(I,J) 

CALL WRTMAT ( DUM »NX , NX»4» • A + B*K «) 

TEMP * 0. 

DO 1 5 I * 1 , NU 
DO 15 J * 1 » NX 

15 TEMP = TEMP + XK( I, J)**2 

TEMO * SORT ( TEMP ) 

WRITFINODEV , 102 ) TEMP 
102 FORMAT (/» 10X , ' I f K J | *»,G15.5) 

C 

C COMPJTF CLOS FD-LDOP EIGFNSOLUT IONS 
C 

I JOB= l 

CALI FIGRF(DUM,NX,4,I JOB , F I G , FVEC , 5, 02 , l ER) 

DO 24 1*1, NX 

FR( I )=REAL( EIG{ I ) ) 

Fid ) = A I MAGCEIGI I )) 

IFIEK I) .LT.O.JGO TO 24 
EMAX = 0. 

DO 20 J * l, NX 

EMAG = CASS (EVECI J, I ) ) 

IF(EMAG.LE.EMAX)GO TO 20 
PMAX = EMAG 

CNORM = CONJG { EVECI J » I ) ) /EMAG* *2 
20 CONTINUE 

DO 22 J = l, NX 

0 1 ( J , X ) * REAL I EVEC ( J , I) *CNORM) 

IF(Enn.GT.O.)Ol(J,dl) = AIMAG(EVECU,I)*CNORM) 
22 CONTINUE 

24 CONTINUE 

DO 30 I * l, NX 
TEMP = 0. 
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00 29 J = l » NX 

28 TEMP » TEMP + Dl(J,n**2 
30 04 ( [ ) = SQRT ( TEMP ) 

C 

C PRINT OUT 6-VALUES AND E-VECTORS 
C 

WRITE(N0DEV,100)NX 

100 FORMAT!/,' THE* ,12,'TH ORDER A+BK .MATRIX HAS El GENSOLUT IONS • , 

1 / , 3X, ' EVAL ( REAL ) • ,4X, 'EVAL! IMAG) ' , 4X , • E- VECTOR * , /) 

DO 32 I = I, NX 

32 WRITEfNODEV, 10 DERI I ) , El ( I) , ( 01 1 J, I ) , J=1 ,NX ) 

101 FORMAT! LX,2G14*6, 10F10.6) 

C 

C COMPUTE AND PRINT OUT THE ANGLES BETWEEN E-VECTORS 
C 

00 38 I = 2, NX 
11=1-1 
00 36 J = 1,11 
T C WP = 0. 

DO 34 K = 1 , NX 

34 TEMP = TEMP + 01 { K , I ) *01 ! K , J ) 

ARG = TEMP/(04(U *D4( J ) ) 

36 03( 1 , J> = 57.2958*AC0S(AMINl(ABS(ARG),Un 

3 3 WRITE(6,104)( I, J,D3! I, J) , J=l, ID 

104 FORMAT!/, 5C ANG* , I 1 , • , ' , 1 1 , ' =*,F7.2)) 

RETURN 

END 

C 

SUBROUTINE XDOT! X, TIME,F ) 

DIMENSION X! 1) ,F(i) 

CO M MU N/ MAT/A! 4,4/ , 0( 4,3 J * XK ( 3 ,4 ) ,NX , NU,NI NT ,NTOT 
C 

C COMPUTE CONTROLS Ull) , U(2)AND 
C X ( 6 ) = INTEGRAL! U( 1 ) **2+U ( 2 ) **2 ) DT 
C F( 6) = U( l)**2+U(2)**2 
C 

TEMP = 0. 

DO 3 I = l , NU 
SUM=0 . 

DO 4 J= l , NX 

4 SUM= SUM+XK { I , J)*X( J) 

XCNINT+I )»SUM 

3 TEMP = TEMP + SUM**2 

F ( N X + 2 ) = TEMP 
C 

C COMPUTE STATES X(L) TO X ! 4 ) 

C AND INTEGRAL OF STATES SQUARED 
C 

TFMP = 0. 

DO 1 1=1, NX 

SUM=0. 

DO 2 J= L » NX 

2 SUM=SUMfA( I ,J)*X( J1 

DO 5 K= 1 , 3 

5 SUM=SUMfB(I ,K)*X{NINT*K) 
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F ( I ) = SUM 

1 TEMP = TEMP + X( I )**2 
F(NX4-1) * TEMP 
RETURN 

END 

C 

SUBROUTINE AUXVAR! X, TIME ) 

DIMENSION X<1) 

C OM MON/ M A T/ A { 4 ,4 ) ,B<4,3) ,XK( 3, A ) , NX , NU.NINT , NTOT 

DO L 1=1, NU 

SUM=0. 

DO 2 J= I , NX 

2 SUM=SUM*XK! I , J)*X( J) 

1 X! N INT+1 )=SUM 

C 

RETURN 

END 

C FILE WRTMAT FORTRAN A1 

C 

C 5/8/79 

C FILE OF UTILITY SUBROUTINES TO SUPPORT VIBE FORTRAN AL 

C 

C WRTMAT - GENERAL MATRIX OUTPUT SUBROUTINE 

C 

SUBROUTINE WRTMAT ( A,N,M, IA ,ANAME) 

0 1 MENS ION A ( I A, l ) 

COMMON/NNOUT/NPRINT 
REAL*8 ANAME 

WRITF(NPRINT, LOO I ANAMF ,N , M 

100 FORMAT!/, ' MATRIX • , A3 , 3X, • ( • , 1 3 , • ROWS X',I3,' COLS)') 
IFfM.LE.LOJGO TO L5 

on 10 i=i , n 

10 WRITF(N PRINT , 101 ) (A( I ,J) , J= 1 , ,M) 

101 FORMAT!/, (IP10E13. 5)) 

RFTURN 

15 DO 2 0 1=1, N 

cl WRITF(NPRINT, 102) C A! I , J ) , J= 1 , M ) 

102 FORMAT! 1P10E13. 5) 

RFTURN 

END 

C 

C TRANSP - TRANSPOSES A MATRIX 

C 

SUBROUTINE TRANSP ! A , N1 , N2 , NA ) 

DIMENSION A ! N A , 1 ) 

COMMON/NNQUT/NOUT 
M = MAXO ( Nl »N2 ) 

IF! M.GT.NA.OR.M.LE.OJGO TO 90 

Ml = M-l 

DO 10 1=1, Ml 

II = I 1 

DO 10 J= I 1 , M 

TEMP = A! I , J) 

A ( I , J ) = A ( J , I ) 

10 A! J, I ) = TEMP 
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RFTURN 

90 WRITE1N0UT, 100)N1,N2,NA 

100 FORMAT ( » *** ERROR IN TRANSP *** Ni,N2,NA = ',JI4) 

RETURN 
END 
C 

C MULT - MATRIX MULTIPLICATION 

C 

SUBROUTINE *ULT ( A , B ,C , N, L , M, NA, NB , NC ) 

DIMENSION A(NA, L) V B(NB,1) ,C(NC, l) 

OOUBLF PRECISION TFMP 
COMMON/NNOUT/NOUT 

IFIN.GT . MINO(NA,NC) . QR.L.GT.NBJGO TO 90 
DO 20 I = l » N 
DO 20 J = 1 , M 
TEMP = 0. 

DO LO K=1,L 

10 TEMP = TFMP ♦ DBLE(AIIfK) > *DBLE t B ( K , J > ) 

20 CII,J) = TEMP 

RFTURN 

90 WRITE INOUT , 100 )M, NA , NC » L» NB 

100 FORMAT!' *** ERROR IN MULT *** N,NA ,NC , L, NB=* , 51 5) 

RETURN 
END 
C 

C M SHI FT - TRANSFEREE. VECTORS AND MATRICES 

C 

SUBROUTINE MSH IFTI A , B , N, M , NA, NB ) 

DIMENSION AINA, 1) , BINS, l) 

COMMON/NNOUT/NOUT 
IF(N.GT.MIN0!NA,N8) ) GO TO 90 
DO 10 1=1, N 
DO 10 J=1,M 

10 B ( I , J ) = A ( I , J ) 

RETURN 

90 WRITE! NOUT , 100 ) N, NA , NB 

100 FORMAT! • *** ERROR IN MSH I FT *** N,NA,NB= • ,31 5) 

RFTURN 
END 
C 

C MATAno - MATRIX ADDITION 

C 

SUBROUTINE M AT ADD! A , B ,C , N , M ,N A , NB , NC ) 

DIMENSION A(NA,L),8(NB,l),C(NC«l) 

COMMON/NNOUT/NOUT 
I F ( N.GT.M INO (NA, NB , NC ) ) GO TO 90 
DO 10 1=1, N 
DO 10 J= 1 , M 

10 C( I,J) = A{ I , J ) * B C I , J » 

RFTURN 

90 WRITE INCUT, 100 ) N, NA , NB , NC 

100 FORMAT! • *** ERROR IN MATADO +** N,NA , MB, NC= ' ,415) 

RETURN 
END 
C 


112 


Appendix C INTODE Computer Program 


C MATSUB - MATRIX SUBTRACTION 

C 

SUBROUTINE MATSUB < A, 8 , C, N , M,NA ,NB ,NC ) 

DIMENSION AINAt 1) ,B(NB, I) tClNCfl) 

COMMON /NNOUT/NOUT 
IF(N.GT.MINQ(NA»N8,NC) ) GO TO 90 
DO 10 1*1, N 
DO 10 J-i.M 

10 C< l,J) = A( 1,0} B ( I » J ) 

RETURN 

90 WRITEINOUT, LOO ) N, NA , NB , NC 

LOO FORMAT ( * *** ERROR IN MATSUB *** N »NA,NB, NC=* , 41 5) 

RETURN 
END 
C 

C SWAP - INTERCHANGES TWO VARABLES 

C 

SUBROUTINE SWAP(A,B) 

C = A 
A = 8 
B = C 
RETURN 
END 
C 

C MS MULT - MATRIX*SCALAR MULTIPLICATION 

C 

SUBROUTINE MSMULT ( S, A , N, M , NA J 
DIMENSION A ( NA , 1 ) 

DO 10 1=1, N 
DO 10 J= l , M 

LO A ( I , J ) = S*A ( l » J ) 

RETURN 

FND 

C 

C ZERO - FILLS A MATRIX WITH ZEROS 

C 

SUBROUTINE Z ERO ( A , N , M , NA ! 

DIMENSION A ( NA, 1 ) 

DO 10 1=1, N 
DO 10 J=l,M 
10 A ( I , J ) = 0. 

RETURN 

END 

C 

c 

c 

C IMAT - LOADS AN ARRAY WITH THE IDENTITY MATRIX 

C 

SUBROUTINE IMAT { A»N,NA ) 

DIMENSION A { N A , 1 ) 

DO 10 1=1, N 
DO 5 J = l , N 
5 A( I , J J = 0. 

10 A < 1,1) = 1 « 

RETURN 


113 


Appendix C INTODE Computer Program 


END 

C 

C 

C ANORM - CALCULATES THE RSS NORM OF A MATRIX 

C 

FUNCTION ANORM ( A* NI , N2,NA) 

OIMFNSION A (NAy I ) 

COMMON/NNOUT/NOUT 
IF ( Ni • GT.NA ) GO TO 90 
ANORM = 0. 

DO LO 1=1, NI 
DO 10 J=l»N2 

10 ANORM * ANORM + A(I,J)**2 

ANORM = SORT { ANORM) 

RETURN 

90 WRITEtNOUT, 100)N,NA 

100 FORMAT ( ' *** ERROR IN ANORM *** M, NA = • ,215) 

RETURN 

END 

C 

C ********** SUBROUTINE PPLOT ********** 

C 

SUBROUTINE PPLOT ( X, Y , NP 0 I NT , XYL I MS , NCHARS ,NL I NE S ) 

C 

C PPLOT GENERATES A PRINTER PLOT UP TO 130 CHARACTERS WIDE 
C BY UP TO 80 LINES HIGH 
C 

C PPLOT INPUTS: 

C 

C XI ) = ARP AY OF X-COORDINATFS OF POINTS 
C Y ( ) = AP RAY OF Y-COORDI NATES OF POINTS 
C NPOINT = TOTAL NUMBER OF POINTS IN CURVES 

C XYLIMSI 4) = LIMITS OF X AND Y AXES IN THE ORDER XMI N, XMAX, YM I N« YMAX. 
C IF XYL I MSI )=0, SCALING IS AUTOMATIC 

C NCHARS = WIDTH OF PRINTER PLOT (NUMBER OF CHARACTERS) 

C NL INFS = LENGTH OF PRINTER PLOT (NUMBER OF LINES) 

C 

DIME NS ION X ( l ) ,Y(1) ,XYLIMS( 1 ) ,XGRID( 12) ,YGR ID( 12) 

INTEGER I FLD( 130,80) , I STAR, IHOR, I VERT, I BLANK, I CROSS 
DATA ISTAR, IHOR, IVERT/'* • , •- ','1 '/ 

DATA I BLANK , ICROSS/ ' *,* + '/ 

COMMON/NNOUT/NODEV 

WRITE (NODEV, 100 ) NPOINT, NCHARS, NL I NES 
100 FORMAT ( / ' PPLOT CALLED ... TOTAL NUMBER OF POINTS =*,Id, 
l //’ NCHARS ,NLINES = • ,214) 

IF( NPOINT. LE.O)RETURN 

IF(NCHARS.LE.130.ANO.NLINES.LE.80)GO TO 1 
WRITE (NODEV, 106) NCHARS.NL INES 

106 FORMAT!/, • LIMITS EXCEEDED IN PPLOT. NCHAR S , NL I NES =*,2110) 

1 XMI N = XYL I MS ( l ) 

XMAX = XYL I MS ( 2 ) 

YM I N = XYLIMSC3) 

YMAX = XYLI*S(4) 

IF(XMIN.GE.XMAX JCALL M INMAX { X » NPQ I NT , XM I N ,XMAX ) 
IF(YMIN.GF.VMAX)CALL MINM AX (Y,N POINT, Y MIN, YMAX) 
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CAU SCALP <XMIN,XMAX,XST,XDEL,NXP> 

XFIN * XST + XDEL*NXP 
OX * XFIN - XST 

1X2 * -XST/DX*NCMARS ♦ l.OOOOOi 

ixaxis * o 

IF! IXZ.GE.l.AND.IXZ.LE.NCHARSHXAXIS * l 
CALL SCALE! YMIN»YMAX» YST , YDEL,NYP) 

YFIN * YST 4- YDEL*NYP 
9Y * YFIN - YST 

IYZ * YFIN/DY*NLINES 4- l.OOOOOl 
I Y AX IS = 0 

IFUYZ.GE.l.AND.IYZ.LE.NLINESHYAXIS * 1 
C WRITE(N0DFV,105)XST,XFIN,YST»YFIN»IXZ,IYZ 

105 FORMAT ( / » • XST, XF IN , YST, YFIN *• t 1P4E 12. 4, 

1 //,* INDEX FOR X ANY Y ZERO REF AXES =*,0P2I6) 

N 1 * NCHARS - 1 
N 2 = NLINFS ~ l 
DO 10 I « 2,Nl 

00 10 J = 2 ,N2 

10 IFLOUyJ) = I BLANK 
DO 20 I = 1, NCHARS 
IF! IYAXIS.EQ.I) IFLD! I, IYZ) ■ IHOR 

1 FLO (1,1) = IHOR 

?0 IFLO( IrNLINES) = IHOR 

DO 30 I = 1,-NLINES 
IF( IXAXIS. EQ.l ) IFLDC IXZ, I) = IVERT 
IFLOUrl) = IVERT 
30 I FLD ( NCHARS , I ) = IVERT 

NX = NXP 4- 1 
NY = NYP + 1 
DO ^5 I = 1 ,NX 
X OP I D( I ) = XST 4- ( I— 1 ) *XOEL 
IX * ( { I-l ) *XDEL/DX)*NCHARS 4- 1.00001 
IX = MINO! IX, NCHARS) 

DO 45 = 1 ,NY 

YGPID(J) = YST 4- < J— 1 ) ♦YDFL 

IY = ( ( J-l ) *YOEL/DY)*NLINES 4- l.OOOOi 

I Y = MI NO ( I Y *NLINES ) 

45 I FLO( IX, IY ) = I CROSS 

WRITE (NODEV, 101 ) (XGRIDI I ) , 1=1, NX) 

L01 FORMAT ( / , * X-GR ID = • , 7F l 0. 3 , / , 9X , 7F 10 .3 ) 

WR I T E ( NQOEV » 102 ) ( YGRIOI I ) ,1=1, NY) 

102 FORMAT!/, • Y-GRID = * , 7F 10. 3 , /, 9X, 7F10 . 3 ) 

00 50 I = 1 , NPO INT 

IX = (!X(I)— XSTJ/DX) * NCHARS + l 
IY = ( ( YF IN-Y ( I ) ) /DY ) *NL INES 4- l 
C I F ( I .LE. 5 0) WRITE (NODEV, 107 ) X ( I ) , Y { I ) , I X, I Y 

107 FORMAT!' X(I),Y!I) =' , IP2E12.4, • I X, I Y =• ,0P2 1 6 ) 
IX = MAXO! MI NO! IX, NCHARS) , l ) 

IY = MAX J !.M INO! IY , NL INES ) , l ) 

5 0 l FLO ! IX, (Y) = I STAR 

nRITF(N0nFV,103) 

103 FORMAT ( l HI ) 

DO 00 IY = 1»NLINES 

60 WRITE (NODFV, 1 04 ) ! I FLD ( I X , I Y ) , I X= l , NCHARS ) 
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104 FORMAT ( IX , 1 30A1 ) 

RETURN 

C DE3UG UNITI6) ,SUBCHK, TRACE, INIT 

C AT 20 

C TRACE ON 

C AT 60 

C TRACE OFF 

END 
C 

C ********* * SUBROUTINE MINMAX ********** 

C 

SUBROUTINE MINMAX ( X, N, XM IN, XMAX ) 

DIMENSION X ( I ) 

XMIN = xm 
XMAX = xm 
IFCN.LE. i) RETURN 
DO 10 1=2, N 

XMAX = AMAX1 ( XMAX, X ( I ) ) 

LO XMIN = AMIN1 ( XM IN, X ( I ) ) 

RETURN 

END 

C 

C********** SUBROUTINE SCALE ********** 

C 

SUBROUTINE SC AL E ( XM I N , XMAX, S TART I , DELI , NPTS I ) 
INTX(X) = IFIX!X-.5*SIGN(*5»X) ) 

I F { XMAX.LE.XMINJGO TO 90 
Xom a XMAX - XMIN 
XEXP = ALOG 10 ( XS I ZE ) 

I C XP = INTX(XEXP) 

XORD = 10.**! EXP 

XNORM = XS I ZE/XORD 

I F ( XNORM • LE . 1 • 6 ) XMOD * .2 

IF ( XNORM. GT • 1 • 6« AND. XNORM. LE. 4. ) XMOD = .5 

I F ( XNORM.GT .4* .AND. XNORM. LE. 8. ) XMOD = L. 

IF! XNORM . GT . 8 . ) XMOD = 2. 

DFLI = XORD*XMOD 
00 10 1=1,30 
XMAG = 10.**( 1-15) 

XPOINT = FLOAT! INTX(XMAG*XMAX)) /XMAG 
IF(X POINT .GE.XMIN)GO TO 20 
10 CONTINUE 

GO TO 90 

20 [SHIFT = (XP0!NT-XMIN)/DEL1 

STAR T 1 = XPOINT - DEL 1*1 SHIFT 

IF( STARTl.GT .XMIN) START l = STARTI - DELI 

NPTS 1 = (XMAX-START1 )/DELl 

IF!STARTl-*-DELL*( FLOAT! NPTSIM-.01 ) . LT . XMAX )NPT S 1 = 
RETURN 

90 rfRITF(6, 100 JXMIN, XMAX 

100 FORMAT!//' ERROR IN SCALE XMIN, XMAX =»,2E12.4) 

RETURN 
END 
C 

c ********** SUBROUTINE APLQT VERSION 1 11/26/80 


NPTS1 + 1 


********* 
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C 

SUBROUTINE APLOTIX, Y, NPT, NPLOT, XYL IMS, WIDTH, HEIGHT, TICKL, 

1 NCASF,XLBL, NXC ,YLBL,NYC,TITL,NTTL, CNAMES, NC 1NDX ) 

C 

C NOTE THP X( ), Y( ) ARRAYS ARF UNCHANGED IN THIS SUBROUTINE 
C 

C APLOT INPUTS: 

C 

C X{ ) *ARR AY OF X-COORDINATFS OF POINTS 
C Y { ) *ARRAY OF Y COORDINATES OF POINTS 
C NPT( )*N0. OF PQ INTS ( X » Y PAIRS) IN EACH CURVE 
C NPL0T=NUMBER OF CURVES 

C XYL IMS( 4)=*LIMITS OF X AND Y AXES IN THE ORDER XMIN, XMAX, YMIN,YMAX 
C IF XYL I MS ( )*0, SCALING IS AUTOMATIC 
C WlDTH=w I DTH OF PLOT IN INCHES 
C HF I 3HT=HEI GHT OF PLOT IN INCHES 

C TICKLE*LENGTH OF TICK MARKS ON AXES IN INCHES AND HEIGHT OF SCALE /VOS 
C IN INCHES. IF TICKL IS NEGET I VF , PUT THE TICK MARKS INSIDE AXES 

C NCAS E*A NO. PRINTED AT THE UPPER RIGHT HAND CORNER OF PLOT FOR 
C INDEXING PURPOSES. 

C XL BL ( ) = ALPHANUMERIC STRING FOR X-AXIS LABEL 

C NXC = NUMBER OF CHARACTERS IN X AXIS LA3EL 
C YLBU ), NYC = LIKEWISE FOR Y AXIS LABEL 
C TITLl ) = ALPHANUMERIC STRING FOR TITLE 
C NTTL » NUMBER OF CHARACTERS IN TITLE 

C CNAMES ( ) = ARRAY OF FOUR-CHARACTER IDENTIFIERS FOR THE NPLOT CURVES 
C NCINDXI ) = INTEGER (.GE.l. AND .LE.NPT ( )) INDICATING LOCATION 
C OF CNAMES IDENTIFIERS FOR EACH CURVE 

C 
C 

DIMENSION X{ L ) ,Y 1 1 ) , XAR( 10 ) ,NPT U ) , XYL I MSI l ) , LSC ALE ( 2 ) 

0 IMF NS I ON XLBL(I) , YLBL ( 1 ) ,T ITL ( I ) , CNAMES! I) , NCINDXI l) 

COMMON/ NNOUT/NODEV 
DATA EPS/I. E-8/ 

WRIT p (NODFV , lOO)NPLOT,NCASE» (NPTI I ), 1=1, NPLOT) 

100 FORMAT!/’ APLOT CALLED NPLOT *• , I 5 , 3X, • NCASE = ' , I 3 , 3X, • NPT ( )»' 
l , 1014, / , ( LOX, 10 1 4 ) ) 
no til) 1 = 1,10 
NNPT = NPT! I ) 

1 10 WR I TE ! NODFV ,109) I , X { I ) , Y ! I +NNPT ) , Y t 1+202 ) , Y ! I +3 03 ) , Y( I +40 4 ) 

109 F0RMATI2X, I 5 , • COORD I NATE S= ' , 6F10.4) 

IF! NPLOT. LE.O)RETURN 
TICK * ABSITICKL) 

HL = 12*T ICK 
H2 = HI + WIDTH 
VI = 9 * T I C K 
V2 = VI + HF I GHT 

FIND TOTAL NUMBER OF POINTS ON PLOT 


NPUINT = 0 

00 2 1=1, NPLOT 

NPOINT = NPQINT + NPT! I ) 

SCALE X-AXES 


117 


n o o *- ooo ooo >- ooo 


Appendix C INTODE Computer Program 


C 

XM IN * XYUMS(l) 

XMAX * XYLIMSI2) 

IF(XYLIMS( l ) .GE.XYLIMS! 2) )CALL MINMAX! X , NPO INT, XM IN , XMAX ) 
CALL SCALE(XMIN,XMAX,XST,XDEL,NXP) 

XFIN « XST f XOEL^NXP 
OX * XFIN - XST 
XF * (H2-H1 )/DX 

WRITE(NODEV,lOl)XMIN,XMAX,XST,XDEL,XFIN,NXP 
101 FORMAT!/* XMINtXMAX =* , IP2E14. 7 , 5X, ' XSTART, XOEL , XFI NI SH 
l 3E 1 0.3 , 5X , • NXP *' ,14) 

SCAIF Y-AXES 

YMIN * XYL IMS ( 3 ) 

YMAX * XYLIMS<4) 

iF(XYLIMSH).GE.XYLIHS(A) )CALL M INMAX ( Y , NPO INT ,YM IN , YMAX ) 
CALL SCALE! YMI N, YMAX, YST, YDEL,NYP) 

YF IN = YST «• YDEL*NYP 
OY = YF IN - YST 
YF = ! V2-V i ) /QY 

WRI TF (NODEV , 102 ) YMI N, YMAX, YST, YDEL , YF IN,NYP 

FORMAT!/* YMI N»YMAX * * , 1 P2E 14 .7 , 5X , • YST ART, YOEL , YFI N I SH *», 
I 3E 10. 3 , 5X , * NYP , 14) 

DRAW OUTER SOROER FOR PLOTS 

CALL PLOT { H l , VI , 3 ) 

CALL PLOT ( H2 , VI , 2 ) 

CALL D LOT ! H2 , V2 » 2 ) 

CALL PLOT !H1,V2,2) 

C ALL PLOT ( H L , VI ,2 ) 

IF APPROPRIATE ADO X = 0, Y = 0 AXES 

IF! XST.GT.O. .OR.XFIN.LT.O. ) GO TO 12 
XX = -XST*XF ♦ HI 
CALL PL0T(XX,V1,3) 

CALL PLOT (XX , V2 » 2 ) 

IF! YST.GT.O. .OR.YFIN.LT.O. ) GO TO 14 
YY = -YST*YF * VI 
CALL PLOT ( HI , YY, 3 ) 

CALL PLOT ( H2 , YY , 2 ) 

CONT INUE 

ADO TITLF, X-AXIS LABEL AND Y-AXIS LA3EL 

HT = TICK 
ANG = 0. 

ANGOO = 

WO - . 9*HT 

= (HlfH2)/2. - • 5 + NTTL * I • 5*w0 
IFINTTL.GT.OJCALL SYMBOL ( XP , V2+4.*HT , l . 5*HT , T I TL , ANG, NTTL ) 
<P = (HL*H2)/2. - .5*NXC*rtD 

IFINXC.ST.O ICALL SYMBOL ( XP , V 1-3 . *HT , HT , XL 3L ,ANG,NXC) 
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YP * <VUV2)/2. - .5*NYC*WD 

I F ( NYC.GT.O JCALL SYHflQL ( Hl-l l. *HT , YP , HT , YLBL » ANG90 * NYC ) 

ADD CASE NUMBER ( IF NONZERO) AND DATE 
NOG * C 

I F INCASE. GT.O JCALL NUMBER (H2*HT »V2»HT » FLOAT I NCASE ) t ANG«NDG) 
CALL 0ATIMP<Hl,V2«-2*HT,HT r »APL0T •) 

ADD X-AXIS TICK MARKS AND SCALE 

N * NXP + L 

YY * VI - 4. *TI CK 

DO 12 1*1, N 

XA * HI 4- ( I- 1 ) *XD6L*XF 
CALL PLOT I XA » VI , 3 ) 

CALL PLOTIXA, V1-TICKL,2) 

XR * XST 4- ( I — l ) *XDEL 
XR 3 XR 4- S IGN { EPS » XR ) 

NDGT * MIND { MAXO l 0*- IF IX ( ALOGIO I ABS I XR ) ) ) +3 ) * B ) 

IF(ABSIXR) .',F. 2.*EPS) N0GT*0 

CALL NUMBER { X A-2 . *T ICK » YY, HT , XR , ANG » NOGT ) 

ADD Y-AXIS TICK MARKS ANO SCALE 

N = NYP 4- 1 

XX = VI - 7 • *T ICK 

DO ? 6 1 3 1 , N 

YA = VI 4- ( I- 1 ) *YDEl*YF 
CALL PL0T(Hl,YA,3) 

CALL PLOTIHl— TICKL»YA»2) 

YY 3 YA 

YR 3 YST 4- ( I — 1 ) *YDEL 
YR 3 YR 4- SIGNIEPS.YR) 

NOGT = MIN0{MAX0(0,-IFIX(AL0GI0(A8SIYR) ) ) 4-3) ,8) 

IFIABSIYR) .LE. 2.*EPS) NDGT=0 
CALL NUMBER (XX,YY,HT»YR,ANG, NDGT ) 

ADD LABEL FOR EACH CURVE 

K = 0 

*0 30 l 3 1 i NPLOT 

I F ( NC I NDX ( I ) «LE .0 ) GO TO 32 

JPT 3 MAXOI MINO(NCINDX( I) ,NPT(I ) ) ,1) 

XX 3 AMINl(AMAXl(Hl,(X(K4-JPT)-XST)*XF4-HI) ,H2) 

YY = AMINKAMAX1 (V1»(Y(K+-JPTJ— YST)*YF4-Vl) » V 2 ) 

CALL PLOT I XX, YY , 3 ) 

CALL °L OT ( XX4-2.*HT , YY4-4. *HT » 2 ) 

CALL SYMBOL (XX4-3. :HT,YY 4-4. *HT # HT,CNAMES( I ),ANG,4) 

N = NP TCI) - I 

RL'T EACH CjRVE 

K 3 K 4- I 

XX = A-UNUAM,AXKHl,{X{K)-*ST)*XF4-HU,H2) 
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30 

C 


10 


20 


90 

100 

c 


10 


YY « AMI NL ( AMAXi ( VI « ( Y( K 
CALL PLOT (XX, YY, 3 ) 

DO 30 J* l »N 
K * K f l 

XX * A M l N l { AM AX L ( H l » ( X { K 
YY * AMINK AMAXUVl, ( Yf K 
CALL PLOT ( XX, YY » 2 ) 
CONTINUE 
R 6 TURN 
END 


l-YSn*YF+VU 


)— XST) *XF *-H I ) 
)— YST) ♦YF+Vl ) 


,V2> 


»H2) 

,V2) 


SUBROUTINE SCALE ( XM IN, XMAX, S TART l , DEL 1 , NP TS l ) 

INTXIX) * IFIXIX— . 5 +SIGNI . 5 ,X) ) 

I F ( XMAX.LR.XMINIGO TO 90 
XSIZE * XMAX - XMIN 
XEXP * ALOG 10 { XSIZE) 

I EXP * INTXIXEXP) 

XORO * 10 ***I£XP 

XNORM * XS I ZE/XORD 

I FI XNORM. LE. 1 . 6 ) XMOD * .2 

I F { XNORM . GT . I . 6. AND . XNORM. LE. 4. ) XMOD - .5 
l F ( XNORM.GT .4. . AND. XNORM.LE.B. )XMOO “ I. 

IF(XNQRM.GT. 8 . ) XMOD * 2 . 

OFLL = XORD*XMOD 

DO LO 1=1-30 

XMAG = 10 .**< I — 15 ) 

X D 0 I N T * FLOAT ( INTX( XMAG*XMAX) ) /XMAG 
T F C XPOINT.GE.XMINJGO TO 20 
CONTINUE 
GO TO 90 

l SHI FT » (XPQINT— XMIN)/DELL 
ST ART l = XPQINT - D£L1*ISHIFT 
IF( STARTI.GT.XMIN)START1 = START 1 - DEL 1 
NPTSI = ( X MAX-START I ) /DEL I 

[F(STARTl+DELl*(FLOAT(NPTSl)+.OL).LT.XMAX)NPTSl = NPTSI + L 

return 

WP. ITEt 6 , 100 ) XMIN, XMAX 

FORMAT!//' ERROR IN SCALE XMIN, XMAX =', 2 E 12 . 4 ) 

RETURN 

END 


SUBROUTINE MI NMAX ( X, N , XM IN , XMAX ) 

o imfnsion xm 

XMIN = X ( l ) 

XMAX = X( 1} 

IF(N.LF.l) RETURN 
DO 10 1=2, N 

XMAX = A MAX I ( XMAX, XI I ) ) 

XMIN = AMINK XM I N , X ( I ) ) 

RETURN 

END 


********** 


SUBROUTINE DATIME 


********** 


SUBROUTINE DATIMEt IMO, l DAY , IYR, I HOURS , I M I N, A.MPM , PNAME , NODEV ) 
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REAi*3 PNAME 

DATA AM / * AM*/, PM/* PM'/ 

CALL 0ATE( IMO, IDAY, IYR) 

CALL STIME( IT IME) 

XHOURS * FLOATt ITIME)/10000. 

AM PM * AM 

IF ( X HOURS. GE. 12. ) AMPM * PM 
I F ( XHOURS .GE • 13 • ) XHOURS * XHOURS - 12. 

I HOURS * XHOURS 

5 XMIN * (XHOURS - IHOURS)*60. 

I M IN * XMIN 

!F( NUDEV. GT .0 ) WRI TE ( NODFV, 100 ) PNAMF, IMO, IDAY , IYR , I HOURS , IMI N , AMPM 
LOO FORMAT ( /' TIME IN ',A8,' IS • ,A2, '/' , A2, • /• » A2, 5X, 12, ' : * , 
l I 2,3X , A4 ) 

RETURN 
END 

***** SUBROUTINE OATIMP ********** 

SUBROUTINE DAT IMP ( XP, YP ,HT, PNAME ) 

RFAL*8 PNAME 
Cn MM ON/NNOUT /NODFV 
ANG * 0. 

WO * .9*HT 

CALL DAT IME ( IMO, IDAY, I YR , l HOURS » l MI N, AMPM ,PNAME » NODE V) 

XHOURS * l HOURS 
XMIN * I M IN 

CALL SYMBOL ( XP, YP »HT, IMO, ANG* 2 ) 

CALL SYMB0L(XP+2.*WD,YP»HT, LH/ , ANG, 1 ) 

CALL SYMB0L(XP+3.*WD,YP,HT, IDAY, ANG, 2) 

CALL SYMBOL ( XP+5 . *WD, YP , HT, LH/,ANG,l) 

CALL SYMBOL ( XP+6.*WD, YP ,HT» IYR, ANG, 2) 

CALL NUMBER (XP+12.*WD,YP,HT, XHOURS, ANG, -L) 

CALL SYMB0L(XP+14.*WD, YP,HT, LH: , ANG, l ) 

CALL NUMBER ( X P+1 5. *WD, YP , HT , XMI N, ANG ,-l ) 

CALL SYMBOL (X P+1 8.*WD,YP,HT, AMPM, ANG, 4) 

RFT'JRN 
END 


FILE INT0D61 INPUT 

INPUT FOR MONTGOMERY'S AIRCRAFT LATERAL SIMULATION 
JUL 2D, 19 80 

£ INTODE 

XI = 0. , . I ,n. , . 1 ,6*0. , 

TMAX “ 5. , 

I PRINT = 2, 

NODFV = 6, 

OT = .325, 

TPP INT = .05 , 

NXPR = l , 2 , 3,4, 5,6,7, 8,2*0, 

A=-3.67994E-1, L.EO ,~2. 4209E-2,2. 53819E-1 ,3*0. ,1.7335E-2, 

“3. 22 79E-2 » 2 • 67949E— l ,-L • 10395E— l »— 9 .6 5926E-* 1, 2 . 61 975E l , 3. , 4. 4ol072E—2 
R=- 7. 6 7 18 3,0., 1.9 6959,0. ,2. 06549, 3., -2. 3 3 943E0, 3. , 
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XK * 1 «‘t941 162E—Q 1,-4. 1923Q29E-0L , 1.0535228E-01 ,2.Q296771E~02» L.3E+00, 
3.6583920E+00,-2. 3151433E-01,-4.7954880E-02, 

NX*4, 

NU*2 * 

(PLOT * 2, 

NXPL » 2, 4, 7, 3, 6*0, 

NCHARS * 80, NUNES » 48, 

NCHARS » 130, NUNES * 62, 

XN4MES*4HR0RA»4HBANK, 4HYAWR , 4H YAW, 4HINX2 » 4HINU2, 4HA I LE » 4HRUDU ,42*0. s 


XLBLSt 1, 1) 

X 

4HR0LL , 

4H RAT, 

4HE (R, 

4HAD/S, 

4HEC ) , 

XLBLS ( L , 2 ) 

X 

4HBANK, 

4H ANG, 

4HLE (, 

4HRAD) , 


XLBL5 (1,3) 

3 

4HYAW , 

4HRATE , 

4H (RA, 

4HD/SE, 

4HC ) , 

XLBLS1 1,4) 

X 

4HYAW , 

4HANGL , 

4HE (R, 

4HAD ) , 


XL3LS1 1,5) 

X 

4HINTE, 

4HGRAL , 

4HXT*X , 



XLBLS (1,6) 

*35 

4HINTE, 

4HGRAL, 

4h(JT*U , 



XLBLSC 1,7) 


4HAILE, 

4HR0N , 

4H (RAD , 

4H ) , 


XL3LSI 1,0) 

** 

4HRUQ0 , 

4HER (, 

4HRAD) , 




NINT * 6, 

NTOT * 3, 

XPLBL * 4HTIME, 4H (SE, 4HC ) ,7*0., 

NXC * 10, 

YPt.BL * 4HSTAT, 4HE , 8*0., 

NYC = 5, 

T[TL*4HM0NT,4HG0ME,4HRY ,4HA/C , 4HLATE , 4HRA L , 4H0YNA ,4HMICS, ,4H K ) 
4H TAU, 4H= , 

TITt 1 13 ) *4H 1.0 , 

NTTL* 52 , 

WIDTH * 5.5, HFIGHT * 2., T I CKL * .06, 

ICONT » 1, 

SEND 
£ INTDDE 

T I TL C 13 ) =4H0 ■ 5 , 

XK = l. 7022794F-01 ,~6 . 01 64034E-01 ,2 .6 19605 l E-0 1 ,-2 . 7018290E-02 , 2 . 9E+00, 
4.053 7797E«-00,-1.1356604E+00,4.89250ldE-01 , 

SEND 

UNTOOE 

T I TL ( 1 3 ) = 4H0 .2 5 , 

XK = 1.23 3632 3E+00*-3. 196 7741 E-Ol ,8. 5040 188E-0 1 , - 1 . 563 7993 E+0 1 , 1 . 9E+00 , 
4.4778833F+00,3.2358761B+0Q»— 1.5113544E+00, 

GEND 
£ 1NTH0E 

TITLl 9)=4H(RC3, 

TITLl 10 ) *4HUST ) , 

T I T L ( 13 )=4H1.0 , 

XK=.6 5 13?, -.335575,. 082589,-. 192543 , 1 . 602 8 5, 3 . 46 768 ,-.73 96,-. 

GFNO 

GINT0D6 

T I TL ( l 3 ) =4H0. 5 , 

XK=. 64359,-. 33107, .6764,-. 244,2. 07949 , 3 . 5064 ,-. 74 1 ,-. 432244 , 

SEND 

UNHIDE 

TITLl 13 ) = 4H0 . 25 , 

XK= 1.0 368, -.3 599, 1.27027, -.322, 2. 05633, 3. 3 1589, -.7* 3469, -.520 

SEND 

SINTC9E 
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TITL(9)*4H(LQR, 

T ITL ( 10)*4H) , RH» 

TITL(li)«4H0* , 

TITU 12 )*4H100. « 

NTTL*48, 

XK*. 3306, -.1784,. 109 18, -.039 l 3, -.69998,. 48 23 7, 2. 6222, -1.4068, 

GENO 
G I NTQ9E 

riTU12)*4Hl. , 

XK» 1.2588,-. 17087, 1.0046,. 0953, -.30636, 1.0 139, 5. 9393, -1.734, 

GFNO 

GINTODE 

ICONT-Q, 

DT* .0 1 , 

TITLI 12)»4H0.Q4, 

XK*5. 17, -.0593 75, 4. 87 19, 1.1309,. 106, 5. 06 l 4,7.2623,-2.5309, 

GENO 

FILE l NT ODE INPUT 

INPUT FOR HALL'S AIRCRAFT LATERAL SIMULATION 
JUL 31, 1980 

GINTODE 

XI * 0 . , . I , 0. , • 1 , 6*0. , 

TMAX * 5, , 

I PRINT * 2, 

NODEV * 6, 

OT = .025, 

TPRINT = .1, 

N 'PR * 1,2, 3, 4, 5,6,7, 8,9,0, 

A=-3. 18, 1. ,-.06, .022,3*0. ,.0644,-63,0. ,-. 2 7,-. 998 ,- 10. 6 , 0. ,4. 5 , 

8*- 14, 4, 3*0. ,1.5,0.,— 2. 59, .037, 2*0. ,-.96,0., 

XK =- 3.9732631 E~ 02,-4.359 8451 E- 03,1. 9987654E— 02 , l .2020696E-02 ,-959E-01, 
2. 5593357F+-00, 2. 16 L 1822E-0 1 ,-3.66 1 126 IE-03 ,4. 75162 79E+00 ,-9. 4E-0 1 , 
-5. 47779200-01,-5. 14 10699E-01, 

NX -4, 

NU=3, 

IPLOT = 2, 

NXPL * 2, 4, 7, 8, 9, 5*0 
NCHARS = 80, NLINES * 43, 

NCHARS = 130, NLINES » 62, 

XNAMES = 4HR0RA , 4HBANK, 4HYAWR , 4H YAW, 4HINX2 ,4HINU2, 4HA I L E , 4HRYAWC , 0 , 


XLSLSt l, 1 ) 

SB 

4HR0LL, 

4H RAT, 

4HE (R, 

4HAD/S ; 

4H EC ) , 

XLBLS ( 1,2) 

S 

4HBANK, 

4H ANG, 

4BL E ( , 

4HRAD ) , 


XL3LS (1,3) 

ss 

4HYAW , 

4HRATE , 

4H (RA, 

4H0/SE, 

4HC ) , 

XL3L S ( 1.4) 

= 

4HYAW , 

4HANGL , 

4HE (R, 

4HAD ) , 


XLBLS ( 1,5) 

*=: 

4HINTE, 

4HGRAL , 

4HXT*X, 



XL3LS (1,6) 


4HINTE, 

4HGRAL, 

4MUT *U , 



XL3L S ( 1,7) 

= 

4HAILE, 

4HR0N , 

4HIRA0, 

4H) , 


XLBLS ( 1,3) 

ss 

4HRUDD , 

4HER ( , 

4HRAD) , 



XLBLS ( 1,9 ) 

= 

4HYAW , 

4HC0NT, 

4HR0L , 

4H ( R AO , 

4H ) , 


NTOT = 9, 

NIMT = 6, 

XPL3L = 4HTIME, 4H (SE, 4HC ) ,7*0., 
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NXC * 10, 

YPLRL * 4HSTAT, 4HE , 8*0., 

NYC * 5, 

TITL*4HHALL,4H AIR»4HCRAF»4HT LA, 4HTERA , 4HL DY,4HNAMI , 4HCSIM, K, 

4H ) T A , 4HU* , 

TI TL ( 12 ) = 4H 1.0, 

NTTL * 48, 

WIDTH = 5.8, HEIGHT * 2., TICKL = .06, 

ICONT « 1, 
aFND 
G INTODE 

TITLI 12)=4H0. 5 , 

XK=. 02469, 00 L394,-. 0 LL66 , . 24055,. 0384, -.0438, -.30 165,. 02299, 4. 14655, 
•OL 14 765, -.3113, 

SEND 

GINTOQE 

T I TL ( 1 2 ) a 4H0. 25, 

XK = 3. 569 5368 E-01 ,-3.949708 l E-02 ,~l .27873 79E-02 , l.4500313E+00,2.9E-0l, 
-8.7618 792E-0l,-l.2922382E*00,2.6900148E-}-00,3.6l001 llE+00,-4. IE+00, 

-1. 245 362 0E + 01, I. 4L50048E+01 , 

SEND 
& INTODE 

TITLI 3)=4HCS(R, 

T I TL( 9)=4H0BUS» 

TI TL ( 10 ) *4HT ) T, 

TITLI I I )=4HAU= , 

TITLI 12J-4H 1.0, 

XK=-. 03489,. L5 L9, -2. 633,. 137833, -.933, 2. 535,. 2 L 75,. 002644, 5. 530 
-.9673 3, -.6965, -.525889, 

SEND 

GINT0DE 

TITLI 1 2 ) = 4H 0.5, 

XK=. 03819, -.00 16245, -.012 548,. 41 33,. 0386 788, -.0443 5, -.27592, 

.0253255, 4. 15544, -.12295, .007238, -.3055, 

SEND 

GINTODE 

TITLI l 2 ) =4H0 .25, 

XK=. 4 9974,. 309 13 l, -3. 344, 1.46229, -.932 129, -5. 79063,. 11 238 7, -.26, 

5.9 79 18, -1.04167 ,—2.05667 ,—1.36633, 

SEND 
G INTODE 

TITLI 8 ) = 4HCS ( L , 

T I T L I 9 ) =4HQR ) , 

TITLI 10)=4HRH0=, 

TI TL I 1 1 ) = 4H 1 00 . , 

NTTL=44, 

XK=.045l 3,. 005 74,. 003798,. 09 16,. 0222, .011679, .056979, .09242, 

.0363 3, -.0754, -.0208, -.008925, 

G P ND 
G INTODE 

TITLI l 1 ) =4H 1.0, 

DT = 0 . 0 1 , 

T 3 LGT=0 . 0 2 , 

XK = . 3653, -.06 7577,. 007833,. 99 586, -.058595, .016246, .1175,. 939 89, 

. 35 148,-. 556 19,-. 27992,-. 1 , 
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SEND 
SINT ODE 

ICONT»0» 

H TL ( II > k 4H0.04, 

£ J"m^r;^:^:' 4 - ,8 *-- 50l5 '-o^’ 7 ..A l eo,,,.r4a. 


-7639, 
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FILE PFRT3 FORTRAN A1 

PROGRAM TO TEST ROBUSTNESS OF LINEAR CONTROL SYSTEMS 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


999 


C 

c 

c 

100 

c 


101 

c 

c 


DOUBLE PRECISION DS6ED 

DIMENSION FVRK4) ,EVI 1(4) ,EVECl(4,4) ,EVR2(4) , EVI2(4) , EVEC2(4,4) 
DIMENSION XX (4004) ,YY(4004) , XYL ( 4 ) , NPTS ( 1 ) 

DIMENSION CNAMES ( 1 ) » XLBL ( LO ) , YL BL ( 10 ) » T I TL( 20 ) ,NCINDX( 1) 

COMMON /MAT/ A ( 4,4) ,B( 4,3 ) * XK( 3, 4) ,CA ( 4,4) ,PA (4,4) , PB ( 4, 3 ) ,N ,M 

COMMON/ I SEED/DSEEDt RAND (28), P 

COMMON/NNOUT/NOUT 

NAMFLIST/PFRT/N,M, ICONT, A,B,XK, P,NSAMP,IPRINT,IPLOT , 

1XLSL , NXC , YLBL ,NYC,TITL, NTTL ,CNAMES,NCINDX 

INPUT VARIABLES IN THE NAMELIST 

N = NO. OF STATE VARIABLES 
M - NO. OF CONTROL VARIABLES 

ICONT = INDEX TO CONTINUE OR STOP MULTIPLE CASE RUNS 
0 , STOP AFTER THIS CASE 
L , CONTINUE TO NEXT CASE 
A = N BY N SYSTEM MATRIX 
B = N BY M CONTROL MATRIX 
XK = M BY N FEEDBACK MATRIX 

P = FRACTION GIVING THF MAXIMUM PFRTUR8ATI0N 
NSAMP = NO. OF SAMPLES 

IPRINT * INDEX TO DETERMINE THE AMOUNT OF PRINTOUT 
I PLOT = INDEX THAT CONTROLS THE GENERATION OF PLOTS 
XLBL = ARRAY CONTAINING THE X-AXIS LABEL FOR PLOTS 
NXC = NO. OF CHARACTERS (LETTERS) IN X-AXIS LABEL 
YLBL = ARRAY CONTAINING THE Y-AXIS LABEL FOR PLOTS 
NYC - NO. OF CHARACTERS IN Y-AXIS LABEL 
TITL = ARRAY CONTAINING A LABEL TO APPEAR IN PLOTS 
NTTL = NO. OF CHARACTERS IN TITLE 

DATA NCASE/O/, ISTART/1/ 

DATA XYL/-10. ,0. ,-3. ,3./ 

DATA WIDTH,HEIGHT,TICKL/ t >.5,3.3,.06/ 

READ (7, PERT) 

DSE6D=1BD0 
NC AS E = NC ASF+ L 
N0UT=6 

NPTS( 1 ) =NSAMP*N+N 
MCUR V= l 

PRINT BANNER 

HRITF(6,100) NCASE 

FORMAT ( 1HL , • ♦♦♦♦PROGRAM PERTURB, CASE 1 , I 5, •****• ) 
WPITE(A.PERT) 

CALL QATI^E( IM, ID, I Y,I H, IMN,AP , 'PERTB *,6) 

W R I T F { 6 , LOl ) TITL 
FORMAT ( / 1 X , 20A4/ /) 

PRINT MATRICES A,9,K AND A +8*K 
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c 

CALL WRTMAT(A,N,N,N,'A* •) 

CALL WRTMAT ( 3,N»M, N, * 3* •) 

CALL WRTMAT<XK,M,N,M,*K= ») 

CALL MATM(A,3,XK,CA, N,,M) 

1 = 1 

CALL WRTMAT(CA»N,N»N» ' A+8*K •) 

C 

C CALCULATE EIGENVALUES ANO EIGENVECTORS 

C 

CALL GENE (C A, N,6VR1, EVIWEVFCl, 1, 1 PR I NT) 
l P = 1 

IF ( I PLOT .GT. 0) CALL XYSTOI XX , YY , EVR1 , EV I 1 ,N, IP ) 

C 

C START MONTE CARLO RUNS 

C 

DO 10 I * L » NSAMP 
CALL PERTB 

CALL MATM(?A,PB, XK,CA,N,M) 

IF( ( I-L)/IPRINT*IPRINT .NE. I-L) GO TO 20 
CALL WRTMAT(PA,N,N,N»'PERT A= •) 

CALL WRT.MAT (P3,N»M,N» 'PERT 3= •) 

CALL WRT.MAT { CA,N,N, N, • PRT A«-BK* ) 

20 CALL OENEICA, N, EVR2» EV12, EVEC2, I » I PR I NT) 

I F { I PLOT .GT. 0) CALL XY STOl XX , YY , EVR2 , EV 12 , N , IP ) 

C 

C CALCULATE ANO PRINT STATISTICS ON REMAX , REMIN ANO RTOMAX 

C 

CALL STATS ( EVRL »EV 1 1 » EVR2 » EV I2t N » I *NSAMP ♦ I PR I NT ) 

10 CONTINUE 

I F ( I PLOT .LE. 1) GO TO 99 
I F ( I START .GE. 1) CALL PLOTSIO, 0, 50 ) 

I F ( I START .LE. 0) CALL PLOT { WI DTH+2 . , 0. ,-3) 

I ST AR T=0 
C 

C 00 PLOTTING 

C 

C ALL APLOT(XX,YY»NPTS?NCURV»XYL,WIOTH,HEIGHT,TICKL,NCASE, 

L XLRL » NXC # YLBL * NYC? TI TL» NTTL » CNAMES t NC I NOX ) 

C 

C GO TO NEXT CASE IF ICCNT IS GRATER THAN ZERO 

C 

99 I F ( I CON T .GE. 1) GO TO 999 

IFUPLOT .GT. 0) CALL PLOT < 0 ., 0. , 999 ) 

STOP 

END 

C 

C MATM - COMPUTES CA = A +■ B * XK 

C WHERE A , 3 AND XK ARE GIVEN MATRICES 

C 

SUBROUTINE MATM( A, 3, XK t CA,N , M) 

0 1 MF NS 1 0 N A ( N » N ) ,3(N»M) , X KIM, N)» CAIN, N) 

00 10 1*1, N 
n n to J- 1 * N 
TEMP=0. 
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no 5 K= l » M 

5 TEMP*TEMP*B{I»K)*XK(K,Jl 

10 CAt t,J) = TEMP«-AU,J) 

RETURN 
END 


PERTB -RANDOMLY PERTURBS THE ELEMENTS OF MATRICES A AND B 

SUBROUTINE PERTB 
DOUBLE PRECISION OSEED 

COMMON/MAT/ A(4 f 4) , B< 4» 3 ) , XK (3 , 4 ) »CA ( 4 ,4 ) , PA ( 4, 4 ) » PB ( 4, 3 ) ,N, M 
COMMON/ I SFE D/DSEED,RAND{ 2BI t P 
NN=N* ( N+M ) 

CALL GOUBSI DSEED»NN*RAND ) 

DO 6 1=1, NN 

RAND ( I )= ( -1 +2*RAND( I ) ) 

NNN=0 

DO l 1=1, N 

PERTURB MATRIX A 

DO 2 J= 1 » N 

NNN=NNN+l 

P A { I , J ) = A C I ,J) 

I F ( A ( I , J ) .EQ. 0. .OR. A ( I , J ) .EQ. L.) GO TO ? 

PA< I ,J)=A( I, J)*( 1. 4-RANO ( NNN) *P J 
CONTINUE 

PERTURB MATRIX B 

DO 3 K=L,M 

NNN=NNN+1 

PtH I , K ) = B ( I , K ) 

I F ( R ( I , K ) .FQ. 0. .OR. B(I,K) .EQ. 1.) GO TO 3 
PB( I * K ) — B ( I « K ) *( 1. f RAND ( NNN )*P) 

CONTINUE 

CONTINUE 

RETURN 

FNO 

STATS - COMPUTES AND PRINTS STATISTICS ON REMAX , REMIN 
AND RTOMAX 

SUBROUT INF STATS < EVRi , EVI l , EVR2* EV I2,N, IS.NSAMP , I PRINT) 

R E AL ♦ B A N AM F 1 , AN AM'S 2 , AN AM E 3 

T IMF NS I ON EVRKN) , EV 1 1 ( N ) ,EVR2(N) , EV 1 2 ( N ) ,R T ( 4 ) 

DIMENSION CREXI 1000) ,CREM( 1000) ,CRTX( 1000), NY 1 20) 

DIMENSION CREXAI 1000 ) , CREMA ( 1000) ,CRTXA( 1000) 

DATA RMAXMA,RMINMA,RTMXMA,SRXA, SRMA , SRT XA /6*0. / 

DATA RMAXM,RMINM,RTMXM,SRX,SRM, SRTX/6*0./ 

DATA ANAMEL , ANAME2, ANAME3/ ' REMAX REMIN RTQMAX •/ 

I F { IS .NF. I ) GO TO 3 

CALL MINMAXIEVRI ,N. REM IN 1 , REMAX 1 ) 

JO l 1=1, N 

ATI I ) =EV 1 1 ( I ) /FVR 1 ( I ) 
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RTMAXL=AMAXl(RT( 1) ,RT ( 2 ) , RT ( 3 ) , RT< 4 ) ) 

3 CALL MINMAX(EVR2,N,REMIN2,REMAX2) 

DO 2 1*1, N 

2 RT( I)=EV12{ n/EVR2(I ) 

RTMAX2*AMAX I(RT( l) »RT(2),RT(3),RT(4) ) 

CREX( I S ) — ( REMAX2-REMAX1 ) /REMAX 1*100. 

CRFM ( I S ) = ( REMIN2-REMI Nl ) /REM INI* 100. 

100 F0RMAT(2X,' PERCENTAGE CHANGES INV/2X, ' REMA X* ' , FI 0-4, //2X, • 
1REMIN=',FL0.4,//2X, 'RTOMAX* ' , F 10,4 ) 

123 F0RMAT(/10X, 'PERCENTAGE CHANGE IN « , A8»/10X , • MAX »,F10.4,' 
IF10.4,' MEAN ',F10.4,' STD DVN *,F10.4) 

CREXA ( I S ) =REMAX2 
CRFMA( IS)=RFMIN2 

124 F0RMAT(2X, 'UNPERTURBED VALUES ARE ' //2X, • REMAX* ', F 10,4, //2X, » 

IREMIN=» ,F10.4,//2X, 'RTOMAX=' , FI 0.4) 

12 5 F0RMATI/10X, • STATISTICS ON ',A8,'WITH PERTURB AT IONS' //10X, ' MAX 

l* »F L0.4» ' MIN', FLO. 4,' MEAN *,F10.4,' STD DVN ',F10.4) 

I F ( RTMAX l .EQ. 0. .AND. RTMAX2 .ED. 0.) GO TO 8 
IF(RTMAX1 .EQ. 0) CRTXI I S ) = 100. 

IF ( R TMAXl .NE. 0.) GO TO 7 
GO TO 9 

7 CRTX( I S 1= ( RTMAX2— RTMAX 1 )/RTMAXl*lOO. 

CRTXA ( IS ) =RTMAX2 
GO TO 9 
C RTX ( 1 S ) — 0 • 

IF ( ( IS-i )/IPR INT*IPRINT .NE. IS-1J GO TO 6 
WRITE (6, 99) REMAXL » REMIN l » RTMAX 1 ,REMAX2, REM IN2,RTMAX2 
99 FORMAT (2X,'RFMAX1 ,REMINL, RTMAX 1 , REMAX2, REMI N2, RTMAX2 =',6E14.6) 
WRITE (6, 100) CREXI IS) ,CREM( IS),CRTX( IS) 

6 l F ( I S .LT. NSAMP ) RETURN 

DO 4 1*1, NSAMP 
RMAXM=RMAXM+CREX< I ) /NSAMP 
R M l N M= R M I NM *CR EM ( I )/NSAMP 
RMAXMA=RMAXMA+CREXA< I ) /NSAMP 
RMINMA=RMINMA+CREMA( I ) /NSAMP 
RTMXMA=RTMXMA+CRTXA ( I ) /NSAMP 

4 RTMXM=RTMXM+CRTX( I ) /NSAMP 
HP 5 1=1, NSAMP 

$RX=SRXMCREX( l)-RMAXM)**2/(NSAMP-l) 

SRM= SRM+ ( CREM ( I J-RMINM) *♦2/ ( NSAMP- 1 ) 

SRXA=SRXA+(CRFXA< I J-RMAXMA) **2/ < NSAMP- 1 ) 

SRMA=SRMA+ ( CREMA ( I )-RM INMA ) **2/ ( NSAMP- l ) 

SRTX A=SRTXA* ( CRTXA ( I )-RTMXMA )**2/ ( NSAMP- 1 ) 

5 SRTX=SRTX+(CRTX( I J-RTMXM ) ♦♦2/ ( NSAMP- l ) 

SRX=SQRT( SRX ) 

SRM=SQRT( SRM) 

SRTX=SQRT{ SRTX) 

SRXA=SQRT(SRXA) 

S R M A= SO R T ( S R M A ) 

SR TXA=STRT( SRTX A) 

CALL M INM AX (CREX, NSAMP, CREXN,CREXX) 

CALL M [NM AX ( CREM,. NSAMP, CREMN,CREMX) 

CALL M INMAX ( CRTX , NS AMP , CR TXN , CR TXX ) 

CALL M INM AX (CREXA, NSAMP, CREXNA, LREXXA ) 

CALL M INMAX ( CREM A, NSAMP ,C RE MNA , CRE.MXA ) 
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CALL M INMAX (CRTXA,NSA.MP,CRTXNA,CRTXXA) 

WR ITE ( 6 * 120 ) 

120 FORM AT ( //2 5X , *ST ATI ST ICS FROM THE MONTE CARLO SIMULATIONS*) 
WRITE (6, 123) ANAME1 , CREXX, CREXN, RMAXM » SRX 

CALL HIST0(CREX,NSAMP,CREXN,CREXX,D,20,NY,NT) 

WR I TE ( 6 » 12 1 ) ANAME1 * D 

121 FORMAT ( l OX , A8 » ' HISTOGRAM , EACH INTERVAL IS *,F10.4»* UNITS*/) 
WRITE (6, 1 22) ( 1,1 = 1,20), (NYU) ,J=1,20) , NT 

122 F0RMATI5X, 2016, //4X, 2016, /40X, 'TOTAL NUMBER OF SAMPLES* • , 16 ) 
WRITE (6, 123 ) ANAME2, CREMX, CREMN.RMINM, SRM 

CALL HISTO ( CREM ,NSAMP ,CREMN ,CREMX ,D , 20, NY, NT) 

WRI TE ( 6 » 121 ) ANAMF2 ,0 

WRITE (6, 122) ( 1, 1=1,20), (NYU), J= 1,20), NT 
WRITE (6, 123) ANAME3»CRTXX»CRTXN ,RTMXM,SRTX 
CALL HIST0(CRTX,NSAMP,CRTXN,CRTXX,0,20,NY,NT) 

WRITE(6,12L ) ANAME3 ,D 

WR IT6 (6, 122 ) ( 1, 1=1,20) , ( NY ( J ) , J=1 , 20) ,NT 
RMAXM=0. 

WRITE (6, 124) CREXA ( IS),CR£MA C I S ) , CRT XA ( IS) 

WRITE (6, 125 ) ANAME1,CREXXA,CREXNA,RMAXMA, SRXA 
CALL HIST0(CR6XA,NSAMP,CREXNA,CREXXA,0,20,NY,NT) 

WRITF(6,12L) ANAME1 , D 

WRITE (6, 122) ( 1,1*1,20), (NY (J),Js 1,20), NT 
wRITF(6, 125) ANAME2, CREMXA , CP EMNA ,RMINMA, SRMA 
CALL HI STO( CREMA, NSAMP,CREMNA,CREMXA, U, 20,NY,NT ) 

WRITE(6,12L) ANAME2 ,0 

WRITE (6, 122) ( I, 1=1,20) , (NY( J) , J=l,20) ,NT 
WRITE* 6, 125) ANAME3 ,CRTXXA,CRTXNA , RTMXMA , SRTXA 
CALL HIST0(CRTXA,NSAMP,CRTXNA,CRTXXA,0,20,NY,NT) 

WRI TE ( 6, 121 ) ANAMF3 , D 

WRITE (6, 122) ( 1, 1=1,20), (NYU ),J= l, 20), NT 
RMAX.MA=0. 

RMINMA=0. 

RTMX.MA=0. 

SRXA=0. 

SRMA=0. 

SRTXA=0. 

RMAXM=0 . 

RM IN M=0 . 

RTMXM=0. 

SRX=0. 

SRM=0. 

SRTX=0. 

RETURN 
FND 

XYSTO - STORES ARRAYS X AND Y IN ARRAYS XX AND YY 

SUBROUTINE X YSTO ( XX , YY , X , Y, N , I P ) 

DOUBLE PRECISION DSEED 

DIMENSION XX ( l) ,YY{ l) ,X(N) ,Y(N) ,RAN0(4) 

DATA SE p O / 1 3 .00/ 

CALI GGU9S( SEED, N, RAND) 

DO 1 1=1, N 

R AND! I )=~l+2*RAN0( I ) 


\ 
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xxup)*xm 

YY ( I P ) = Y ( I J 

l F ( Y ( I) .EQ. 0.) YY ( I P ) = YY( IP) +RAND ( I ) *• 04 
I IPs-IP^-l 

RFTURN 
END 
C 

C H I STD - COMPUTES HISTOGRAM 

C 

SUBROUTINE H I STO ( X, N, XM I N ,XMAX, D ,NH, NY , NT ) 

DIMENSION X { l) ,NY(NH) » XH( 20) 

D=( XMAX-XMINI/NH 
XDUM=XMI N 
DO l 1=1, NH 
NY ( I )*0 
X H ( I )=0+XDUM 

1 XDUM=XH( I ) 

XH(NH)=XMAX 
DO 2 1=1, N 
DO 3 J= I , NH 

I F ( X ( I) .LF. XH( J) ) GO TO 2 

3 CONTINUE 

C5 I FC J .EQ. I .OR. J .EQ. 20) WRlTE(6,10l) X(I) 

2 NYI J)=NY( JM-I 

101 FORMAT! 15X, 5FI0. 4) 

NT=0 

DO 4 1 = L , NH 

4 NT=NT4-NY(I) 

RETURN 

END 

C FILE GENE FORTRAN 4/23/80 

C 

C PROGRAM TO FIND El GENSOLUT IONS OF A GENERAL MATRIX 
C 

SUBROUTINE GENE ( A, N , EVR, EVI , E VEC , I S , I PRT ) 

D IMF NS I ON A (4,4) ,EVR(4) , EV I ( 4) , EVEC ( 4,4 ) , WORK! 25) 

commqn/nnout/nout 

NOUT = 6 

C CALL WPTMATt A, 4, 4, 4, • A •) 

N = 4 
IA = 4 
MODE = I 
IE = 4 
I PR t NT = 2 

CALL P I GER ( A, N, l A, MODE, EVR, EVI , EVEC, I E, WORK , I PRINT, *A + B*K 
II S, IPRT ) 

RETURN 

END 

C 

SUBROUTINE ESHIFT ( CEVAL , CEV6C ,N ,NC , EVALR , EV AL I , EVEC , WORK ,NE ) 

COMPLEX CEVAL (1) ,CEVEC(NC, l ) ,CNCRM 

DIMENSION EVALR ( l ) , FVALI ( 1) ,EVFC(NE, l) ,W0RK(N6, 1) 

DIMENSION I N 0 X ( 16) , EVALAt 16) 

DO 10 1=1, N 

EVALR ( I ) = REAL ( CEVAL ( I ) ) 
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EVALI(I) = A I MAG( CEVAL ( I ) ) 

EVALA( I ) * CABS (CEVAL ( I ) ) 

IF ( EVAL I ( I ) )26, 15, 15 
15 EMAX * 0. 

DO 20 J= 1 » N 

EMAG * C ABS( CEVEC ( J » I ) ) 

I F C EM AG. LE. EM AX) GO TO 20 
FMAX = FMAG 

CNORM * CONJG(CEVEC(J#I ) )/£MA0**2 
20 CONTINUE 

DO 22 J=1,N 

22 WORK I J » I ) * REAL (CEVEC ( J # l ) *CNORM) 

GO TO 10 

26 DO 28 J= l » N 

23 WORK(J,I) * AIMAG(CEVEC(J,I-l)*CNORM) 

10 CONTINUE 

CALL ORDER ( EVAL A# I NOX# N» 2,£vALR»EVALI ) 

00 30 J= 1 » N 
JJ = INDX(J) 

DO 10 1*1, N 

30 FVEC(IfJ) * WORK ( I # J J ) 

1 = 0 

40 1 = 1 + 1 

! F C r .GT.NLRETURN 
I F ( EVALI ( I ) ) 45,40,45 
45 EVALIU) = ABS ( EVAL 1(1)) 

1=1 + 1 

FVAL MI) = -FVALI ( 1-1 ) 

GO TO 40 
END 
C 

C ORDER - ORDERS THE ELEMENTS OF A VECTOR 3Y ALGEBRAIC SIZE. THE 

C RESULTING PERMUTATION OF THE INDICES IS RETURNED IN VECTOR TV. 

C I*= MODE = 1, THE ELEMFNTS OF VECTOR VI ARE ALSO REORDERED 

C IF MODE = 2, 30TH VECTORS VI AND V2 ARC REORDERED 

C 

SUBROUTINE ORDER l RV , I V , N , MODE , VI , V2 ) 

DIMENSION RV( l) ,IV( U ,V1( l) ,V2( 1) 

I F( N.LE. t ) RETURN 
DO 1 1=1, N 

1 I V ( I ) = I 

DC 4 II = 2 » N 
I = I I 

2 J = I. - 1 

IF (RV( I) .GE.RV( J) )G0 TO 4 
CALL SWAP (R V( I ) ,RV( J) ) 

I p { MODE.GE. I ) CALL SWAP ( V l ( I ) , VI ( J ) ) 

I F ( MODE.GE. 2 ) CALL SWAP ( V2 ( I ) , V2 ( J ) ) 

KK = I V ( I ) 

I VC I ) = IV( J ) 

IV(J) = KK 
IFU.LE.DGC T0 4 
1 = 1-1 
GO TO 2 

4 CONTINUE 
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RETURN 
END 

EIGER - CALCULATES THE EIGENVALUES OF A MATRIX ( MODE * 0) 

OR ROTH EIGENVALUES AND EIGENVECTORS I MODE * l) 

THE EIGENVALUES ARE ORDERED BY ABS VALUE 
[PRINT CONTROLS THE PRINTOUT 

[PRINT * L , LIST EIGENVALUES ONLY 
l PR I NT * 2 f ALSO LIST EIGENVECTORS 

SUBROUTINE E I GER ( A, N, I A, MODE , EVR , EV I , EVEC , I F ■, WORK , 1 PR I NT, ANAttff, 
IIS, IPRT) 

DIMENSION A ( I A, 1 ) , EVECI I E, l ) , EVR ( 1 ) , EV I (1) , WORK ( 1 ) 

REAL*8 ANAME 

COMPLEX CEVAL! 16), CEVEC! 16, 16) 

COMMON/NNOUT/NPRINT 
IF(N.GT.I6)ST0PI 
I JOB * MODE 
IC * 16 

CALL FIGREI A,N, IA, I JOB,CSVAL,CEVEC, I C, WORK, IER) 

I FI I ER.GT.0)WRITE(NPRINT,102) IER 
10? FORMAT!/, • ERROR IN EIGRF IER**, 15) 

IF! IPRINT.GE.3)CALL CEVPRTI CEVAL ,CEVEC ,N, IC »WGRK,MODE) 

CALL ESH I FT ( CEVAL , CEVFC , N , IC, FVR, EV I » EVEC ,wORK, IE) 

I F ( IPRINT.LE.O) RETURN 

[ F ( I PR INT . EO. I ) GO TO 40 

IF! I IS- 1 )/IPRT*IPRT .NE. IS-I) RETURN 

IFIMODF.LF.OIGO TO 40 

WRITEINPRINT, 104) N, ANAME 

104 FORMAT!/' THE', 1.3,' ORDER »,A8,' MATRIX HAS', 

l /,3X» 'EVAL(REAL) ' ,4X, 'EVAL! I MAG)' , 4X , • E-VFCTOR ' ) 

DO 30 1=1, N 

30 WRITE ! NPR INT , 103 ) EVR! I ) , EV I { I ) , !EVEC! J , I ) ,J=l,N) 

103 FORMAT! IX, 1 P2 El 4. 6 , OP l OF 1 0. 6 , / , !29X, 10F10.6) ) 

RETURN 

40 WRITEINPRINT, 105) N, ANAME, ! EVR ( I ) , 1= 1 , N ) 

105 PUPMAT!/' THE', 13,' ORDER *, A8, ' MATRIX HAS EIGENVALUES' 

l ,' (REAL)/! IMAG)' , / ! I P 1 0E1 3 . 5 ) ) 

WRITEINPRINT, 106) (EVI ( I) , I=1,N) 

106 FORMAT! /! IP10E13-5) ) 

RFTURN 

END 

C 

SUBROUT INE CEVPRT(CEVAL,CEVEC,N, I C , WOR K ,MGOE ) 

COMPLEX CEVAL ( l ) , CEVEC ! IC , l ) 

DIMENSION WORK { N, 1 ) 

DO LO 1=1, N 

WORK (1,1) = REAL (CEVAL ! I ) ) 

ID WORK! 2, I) = A I MAG (CEVAL! I ) ) 

CALL WRT MAT ( WORK , 2 »N» N , ' E-VAL S •) 

DO 20 1=1, N 
DO 20 J= 1 , N 

tQ WORK (I, J) = REAL ! CEVEC ( I , J ) ) 

CALL wRT MAT ( wORK».N , N» N» ' RL E-VEC * ) 

DO 30 1=1, N 
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00 30 J*l,N 

30 WORK { I , J ) * AIMAG!CEVEC ! I » J H 

CALL WRTMAT(WORK,N,N,N, » I* E-VEC' ) 

RETURN 

END 

C 

C WRTMAT - WRITES MATRIX A 

C 

SUBROUTINE WRTMAT! A, N*M» IA»ANAME) 

DIMENSION A! IA,i) 

COMMON/NNOUT/NPRINT 
REAL*B ANAME 

WRITEfNPRINT, 100 ) ANAME » N » M 

100 FORMAT!/, • MATRIX • » AB , 3X » M ' » I 3 , ' ROWS X*,I3,' COLS)') 
IF(M.LE.10)G0 TO 15 

00 LO 1=1, N 

10 WRITE! NPRINT,101)(A{I,J),J X 1,M) 

101 FORMAT!/, { IP10E13. 5) ) 

RFTURN 

15 DO 20 1 = 1, N 

20 WRITE INPRINT, 102) !A( I ,J) ,J=l,M) 

102 FORMAT! 1P10E13. 5) 

RETURN 

END 

C 

C M SHI FT - TRANSFERES VECTORS AND MATRICES 

C 

SUBROUTINE MSHIFT { A , B , N»M ,NA , N8 ) 

01 MANSION A ( N A , 1 ) ,B!NB, 1) 

CQMMON/NNOUT/NOUT 
IF!N.GT.MINO(NA,NB) ) GO TO 90 
00 10 1=1, N 
DO 10 J= 1 , M 

10 B( I , J ) = A( 1,4) 

RETURN 

90 WRITE! NOUT , 100 ) N* NA , NB 

100 FORMAT! • *** ERROR IN MSHIFT *** N ,NA ,NB= • , 31 5) 

RETURN 

END 

C 

C SWAP - INTERCHANGES TWO VARABLES 

C 

SUBROUTINE SWAP!A,8) 

C = A 
A = B 
B = C 
RETURN 
END 
C 

C********** SUBROUTINE APLOT ********** 

C 

SUBROUTINE APLOT (X, Y, NPT ,NP LOT, XYLIMStwIDTH, HEIGHT ,TICKL , 
l NCA$f=,XLBL» NXC , YLBL , NYC , TITL ,NTTL , CNAMES , NC INDX ) 

C 

C NOTE THE X! ), Y( ) ARRAYS ARE UNCHANGED IN THIS SUBROUTINE 
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C 

C APLOT INPUTS: 

C 

C X( ) * AP RAY OF X-COORDINATES OF POINTS 
C Y ( ) = ARRAY OF Y COORDINATES OF POINTS 
C NPT ( )*NO. OF POINTS!X,Y PAIRS) IN EACH ClRVE 
C NPLOT *NUMBER OF CURVES 

C XYLIMSI 4)*LIMITS OF X AND Y AXES IN THF ORDER XMIN, XMAX, YMIN»YMAM 
C IF XYLIMSI )*0, SCALING IS AUTOMATIC 
C WIDTH*WIDTH OF PLOT IN INCHES 
C HEIGHT* HEIGHT OF PLOT IN INCHES 

C TICKLE*LENGTH OF TICK MARKS ON AXES IN INCHES AND HFIGHT OF SCALE A-OS 
C IN INCHES. IF TICKL IS NEGET I VE , PUT THE TICK MARKS INSIDE AXES 

C NCASE*A NO. PRINTED AT THE UPPER RIGHT HAND CORNER OF PLOT FOR 
C INDEXING PURPOSFS. 

C XLBLI ) = ALPHANUMERIC STRING FOR X-AXIS LABEL 
C NXC * NUMBER OF CHARACTERS IN X AXIS LABEL 
C YLBLI ), NYC * LIKEWISE FOR Y AXIS LABEL 
C TITLI ) = ALPHANUMERIC STRING FOR TITLE 
C NTTL = NUMBER OF CHARACTERS IN TITLE 

C CNAMES( ) = ARRAY OF FOUR-CHARACTER IDENTIFIERS FOR THE NPLOT CURVES 
C NCINPX! ) » INTEGER UGE.l. AND .LE.NPT! )) INDICATING LOCATION 
C OF CNAMES IDENTIFIERS FOR EACH CURVE 

C 
C 

0 l MENS ION X(l),Y(l) ,XAR( 10) ,NPT!l) .XYLIMSI l ) 

DIMENSION XLBL(l) »YLBL( l) ,TITH l) , CNAMES ( 1) ,NCIN0XU) 

CCMMON/NNQUT /NGOEV 
DATA EPS/l.E-8/ 

WRITEINODEV, l 00 ) NPLOT ,NC ASE , ( NP T ! I ) , 1= l , NPLOT ) 

100 FORMAT!/' APLOT CALLED NPLOT = • , I 5 , 3X , • NCASF » ' , I 3, 3X . • NPT ( ) a' 

1 , 1014,/, I 10X, 1014) ) 

DO 110 1=1,10 
NNPT = NPT (I ) 

110 WRITEINODEV, 109) I,X(I),Y( I *NNPT ) , Y { I + 202) , Y( l *303 ) , Y(I *404) 

109 FORMAT I 2X, 15, 'COORDINATES** « 6F10.4) 

IF! NPLOT. LE.OJRETURN 
XMN = XYLIMS! 1) 

XMX = XYL IMS l 2 ) 

YMN = XYLIMS! 3) 

Y MX * X Y L I M S ( 4 ) 

TICK = ABS(TICKL) 

HL = L2*T ICK 
H2 = HI * WIDTH 
VI = 9*T ICK 
V2 = VI * HEIGHT 

FIND TOTAL NUMBER OF POINTS ON PLOT 

NPOINT = 0 
DO 2 1=1, NPLOT 
NPOINT = NPOINT + NPT! I ) 

IF! XMN.GE.X.MX )G0 TO 4 
XMIN = XMN 
XMAX •= XMX 
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00 TO 6 

4 CALL M INMAX ( X »NPO INT, XMI N»XMAX ) 

6 IF<YMN.G£.YMX)GO TO 8 

YMIN * YMN 
YMAX * YMX 
GO TO 10 

8 CALL MINMAX(Y,NPOINT, YMIN, YMAX) 

C 

C SCALE X-AXES 
C 

10 CALL SCALE(XMIN,XMAX,XST,XOEL,NXP) 

XFIN = XST ♦ XDEL *NXP 
OX * XFIN - XST 
XF - (H2-H1J/DX 
C 

C SCALE Y-AXES 
C 

CALL SCALE (YMIN, YMAX, YST , YDEL ,NYP ) 

YF IN * YST + YD£L*NYP 
OY * YF I N - YST 
YF * {V2-V1J/DY 

WRITEINOOEV, 101 ) XMI N,XMAX,XST, XDEL, XFIN,NXP 
WR I TF ( NODEV , 1 02 )Y MIN, YMAX, YST, YDEL, YFIN ,NYP 

101 FORMAT!/* XMIN,XMAX = * , IP2E14.7 , 5X , * XST4RT, XDEL , XFI N I SH * • , 
1 3Ei0.3,SX,'NXP =• , 14) 

102 FORMAT ( / * YMIN, YMAX * * , 1 P2F 14. 7 ,5X, ' YSTART, YDEL , YFI N I SH = * , 
1 3f 10. 3 , 5X » ' NYP *• , 14) 

C 

C DRAW OUTER 30R0FR FOR PLOTS 
C 

CALL PL0T(HL,V1,3) 

CALL PL0TIH2, VI, 2) 

CALL PL0T(H2,V2,2) 

CALL PLOT (III ,V2 ,2 ) 

CALL PLOT ( H 1 , VI ,2 ) 

C 

C IF APPROPRIATE ADD X = 0, Y * 0 AXES 
C 

77 IF ( XST.GT.O..OR.XFIN.LT.O. ) GO TO 12 

XX = - XST*XF + HI 
CALL PLOT (XX, Vi ,3) 

CALL PLUT(XX,V2,2) 

12 IF(YST.3T.O..OR.YFIN.LT.O.)GO TO 14 

YY = -Y$T*YF Vi 
CALL PLOT ( H l , YY , 3 ) 

CALL PL0TIH2, YY,2) 

14 CONTINUE 
C 

C ADD TITLE, X-AXIS LA3EL AND Y-AXIS LABEL 
C 

HT = TICK 
ANG - 0. 

ANG90 * 90. 
rfO = HT 

XP * (H1+H2J/2. - .5*NTTL*1 .S*WD 


136 


OOON O OOOru O OOO O O O 


Mppenaix u PtRTB Computer Program 


l F (NTTL .GT • 0) CALL SYMBOL ( XP , V2*4. *HT , i . 5*HT , T I TL, ANG, NTTL) 
XO . (Hl+H2)/2. - • 5*NXC*W0 

I F ( NXC. GT . 0 )C ALL SYMBOL ( XP, V 1-8, *HT , HT , XLBL , ANG, NXC ) 

YP * (Vl+V2)/2. - • 5*NYC*WD 

[ F { NYC .GT .0 )C ALL SYMBOL ( HI- l l .*HT , YP ,HT , YLBL , ANG90 ,NYC ) 

A00 CASE NUMBER (IF NONZERO) AND DATE 
NOG * 0 

IF(NCASE.GT.O)CALL NUMBER (H2+HT, V2, HI , FLOAT! NCASE ) , ANG, NOG) 
CALL OAT I MP ( HI »V2+2*HT »HT»'APLOT •) 

A 00 X-AXIS TICK MARKS AND SCALE 


N * NXP + 1 

YY * VI - 4. *TICK 

00 22 I*l,N 

XA * ( I- 1 ) *XOEL*XF 

CALL PL0T(HUXA,V1,3) 

CALL PLOT(HlfXA,VI-TICKLf 2) 

XX * HI XA - 2.*TICK 
XR = XST + ( I- l ) ♦ XDEL 
XR * XR SIGN(EPS,XR) 

IF { ABSI XR) .LE«2v»*EPS)XR = 0. 

NQGT = m(NO(MAXO(G,-IFIX( ALQGLO { ABS { XR ) ) ) +3 ) 
I F ( A3S ( XR ) .LE. 2.*EPS) NDGT=0 
CALL NUMBER (XX, YY,HTtXR, ANG, NDGT) 

CONTINUE 


ADO Y-AXIS TICK MARKS AND SCALE 

N x NYP + L 

XX x VI - 7. ♦TICK 

00 24 1= I , N 

YA = ( I- I ) *YDEL*YF 

CALL PLOT(Hl, Vl+YA,3) 

CALL PLOT( Hk-TICKLtVl+YA, 2 ) 

YY * VI + YA 

YR = YST ♦ ( 1-1 )+YDEL 

YR * YR 4- S IGN(EPS,YR) 

IF( \9S ( YR ) • LE • ♦EPS ) YR = 0. 

NDGT = MINO(MAXO(O,-IFIX(ALOG10(ABS(YR) ) )+3) ,8) 
IF { A BS ( YR ) .LE. 2.*EPS) NDG T=0 
CALL NUMBER (XX, YY,HT,YR, ANG, NDGT) 

CONTINUE 


ADO LABEL FUR EACH CURVE 
K = 0 

DO 30 1=1 , NPL3T 

IFtNCINOXm .LE.01G0 TO 28 

JPT * MAXO( MIN0( NC I NOX( I ) ,NP T ( I ) ) , I ) 

XX = A MINI IAMAX1 (HI, (X(K+JPT)-XST)*XF+H1) ,H2) 
YY = AMI Nl( AMAXl (VI , ( Y( K + JPT l-YST )-*YF + VI ) ,V?) 
CALL PL0T(XX,YY,3) 
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CALL PLOT(XX+2.*HT,YY+4.*HT,2) 

CALL SYMBOL (XX +3 .*HT,YY+4.*HT»HT,C NAMES! I),ANG,4 
28 N * NPT( U - 1 

PLOT EACH CURVE 

K * K + 1 

XX = AMINK AMAX1!H1,(X(K)-XST)*XF+H1) ,H2) 

YY « AMINK AMAXKV1,(Y(K)~YST)*YF+Vl) ,V2) 

CALL PL0T(XX,YY,3) 

C CALL SYMBOL (XX»YY».04»1H.»0.»1) 

CALL PLOT ( XX , YY » 2 ) 

CAL'. CIRCLE(XX,YY,.02) 

00 30 J= l »N 
K = K + 1 

XX = AMINl! AMAXl (H! ,(X(K)-XST)*XF+Hl) ,H2) 

YY = AMINKAMAXl(Vl,(Y(K)-YST)*YF+Vl),V2) 

CALL PLOT(XX,YY, 3) 

CALL PLOT ( XX , YY » 2 ) 

C CALL SYMBOL (XX, YY,.04»1H.»0*,1) 

I F ( J .LE. 3) CALL C I RCLE ( XX , YY , . 02 ) 

30 CONTINUE 

RETURN 
END 
C 

SUBROUTINE SCALE! XMIN»XMAX» START1 »DEL l »NPTS l ) 
INTX(X) = I FI X ( X— • 5+SIGN ( • 5 » X) ) 
IF(XMAX.LE.XMIN)GO TO 90 
XSIZE = XMAX - XMIN 
XEXP = ALOG 10 ( XSIZE) 

■i.:/’’ = INTX(XEXP) 

XOfiC = 1Q.**I EXP 

XN‘>* M * XSI ZF/XORD 

IF(X NORM. LE.1.6)X MOD = .2 

I F ( XNORM. GT • I • 6 .AND • XNOP M. LE • 4. ) XMOD = .5 

IF(XN0RM.GT.4..AND.XN0RM.LE.8.)XM0D = 1. 

IF ( XNORM.GT .8. ) XMOD = 2. 

DELI = XORD*XMOD 
DO 10 1=1,30 
XMAG = L0.**( I — L 5 ) 

XPOINT = FLOAT! I NTX ( X MAG*XMA X ) ) /XMAG 
IF! XPOINT.GE.XMIN) GG TO 20 
10 CONTINUE 

GO TO 90 

20 [SHIFT = ( XPO I NT-XMI N) /DEL I 

START1 = XP Cl NT - 0F:,1*I SHIFT 

1 F ( START1 . GT .X.M IN ) START L = START! - DELI 
NPTS1 = ( XMAX- ST ART 1 ) /DEL 1 

IF ( START l + OELK! FLOAT! NPTS1 ) +.01 ) . LT. XMA X ) NPT S 1 
RETURN 

90 WRITEI6, 100JX MIN, XMAX 

100 FORMAT!//’ FRROR IN SCALE XM I N , XMA X = • , 2E 12 . 4 ) 
RETURN 
END 
C 


•NPTS1 + 1 
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SIJBROUT INF MI NMAX(X»N, XMIN, XMAX) 

DIMENSION X(I) 

XM IN = X(I) 

XMAX = X C 1 ) 

IF(N.LF.l ) RETURN 
DO 10 1 = 2, N 

XMAX = AMAX 1 ( XMAX » X ( I ) ) 

10 XMIN = AM INI { XM IN, X ( I ) ) 

RETURN 

END 

********** SUBROUTINE OATIME ********** 

SUBROUTINE OAT I ME ( IMO, IDAY, IYR, I HOURS, I M l N, A^PM, PNAME ,NOOEV ) 

REAL*9 PNAMF 

DATA AM/' AM'/, PM/' PM'/ 

CALL DATE! IMO, IDAY, IYR) 

CALL STIMF(ITIME) 

XHOURS = FLOAT! ITIME)/10000. 

AMPM = AM 

IFIXH0URS.GE.12. )AMPM = PM 
I F ( X HOURS « GE • 13 • ) XHOURS = XHOURS - 12. 

I HOUR S = XHOURS 

XMIN = (XHOURS - I HOURS ) *60 ■ 

I M I N = XMIN 

IF ( NODEV. GT.O )WRI TE (NODEV , 100) PNAME, I MU, lyAY, IYR, ! HOURS , I M IN , AMPM 
00 FORMAT!/' TIME IN ',A8,» IS • , A2 , ' / ' , A2 , ' /' , A2, 5X , 1 2 , • : • , 

1 I 2 » 3X » A4 ) 

RETURN 

END 

********** SUBROUTINE OATIMP ********** 

SUBROUTINE DAT I MP ( XP , YP , HT , PNAME ) 

REAL *8 PNAME 
COMMON/NNOUT/NODEV 
ANG = 0. 

WD = .9*HT 

CALL DATIME! IMO, IDAY, IYR, IHOURS , IMIN, AMPM, PNAME , NODEV) 

XHOURS = IHOURS 
XMIN = IMIN 

CALL SYMBOL <XP,YP,HT, I MO, ANG, 2) 

CALL SYMBOL ( X P+2. *WD, YP , HT, IH/,ANG,1) 

CALL SYMBOL (X P+3. *WD,YP,HT, I DAY, ANG, 2) 

CALL SYMBOL (XP + 5.*WD, YP,HT, IH/,ANG,1) 

CALL SYMBOL l X P+6. *WD» YP * HT, IYR, ANG, 2) 

CALL NUMBER!XP+12.*WD,YP,HT, XHOURS, ANG,- l ) 

CALL SYMBOL (X P+14. *WD, YP , HT , 1H: , ANG, l I 
CALL NUMBER(XP+15.*WD, YP , HT , XM I N , ANG ,-l ) 

CALL SYMBOL !XP+18.»WD,YP,HT, AMPM, ANG, 4) 

RETURN 

END 

PILE PERT INPUT 

IT CONTAINS THE INPUT FOR THE MONTGOMERY AIRCRAFT 
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c 

£PERT 

N*4, 

M* 2 i 

ICONT=l, 

I PR INT* 1000, 

IP LOT = 2 » 

XLBL* 4HRE AL » 4H( LAM, 4HDA ) , 

YLBL=4HIMAG»4H INAR, 4HY ( l A, 4HMDA) , 

NXC* L 2 » NYC»16, 

TITL=4HM0NT,4HG0ME,4HRY A,4H/C C,4HL0SE,4HD L0»4H0P E, 4HI GEN, 44HE PE, 
4HRTUR , 4H8AT I ,4H0NS, ,4HP=. I ,4H |MI,4HN K),4HTAU=,4H l.Q* 
NTTL=72,CNAMES=4H ,NCIN0X=0, 

NSAMP»5, 

P=. 10, 

A=-3. 679 8 4E-1, l .EO, -2.4209 E- 2,2.53819 E-i, 3*0. ♦ 1.7835 E-2» 

-3. 2279 E- 2,2. 6 79496- 1, — 1.10395 E—l ,-9. 659266- l , 2 .6 l 875EI , 0. ,4.46 107 2E-2 
8= -7. 67183,0., 1.96959,0. ,2.06649,0. ,-2. 33843F0, 0. , 

X K * 1.4941 162E-0 l ,— 4. 1928029E-0 1 , 1 . 0535228E-01 , 2 . 029677 IF-02 , 1. 6E+0 0 , 

3. 65839206+00,-2.31514336—01 ,-4. 79543805-02, 

£FND 

SPERT 

TI TL ( 18 ) =4H0 • 5 , IC0NT=1, 

XK » 1. 70 22 79 4F— 01,-6.01640345-01 , 2. 6 196051 E-01 , -2 . 7 01 3290 E-02 ,2 .9F+00, 

4. 0587797F + 00,— l. I 85o604 E+00 ,4.8925018 F— 0 1 , 

SEND 
CPFR T 

TI TL ( 13 ) = 4H0. 25, 

XK = 1.23 36323E+00, -3.1967741E-0l*8. 5040188 E— 01,-1.563 79 93 E+ 01,1. 9E+ 00, 
4. 47 7 883 3E+0 0,3. 235876 1E + 00,- 1.5 1135 44E +00, 

6FND 

£PERT 

TI TL ( 15 J =4Hi { RQ, T I TL ( 1 6 ) = 4HBUST , T I TL ( 17 )=4H) , TA,T I TL ( 13) = 4HU= 

TITU 19 ) = 4H 1.0, 

NTTL=76, 

XK = . 6 5 13 2, -.3 8 55 75,. 08 25 89,-. 192 548 , 1 . 602 8 5, 3 . 46 76 8 ,-. 7396,-. 60 73T , 

SEND 

SPERT 

TI TL ( 19 ) = 4H0. 5 , 

XK=. to 4359, -.3 310 7,. 6764, -.244, 2. 07949, 3. 5064, -.74 l, -.432 244, 

£ END 
£PERT 

TI TL ( L9)=4H0.25, 

XK=l. 0 3 68, -.3 599, 1.2702 7, -.322, 2. 0568 3, 3. 3 1539,-. 74 3469 520 8 6, 

£ END 
£P6PT 

TITU 15 ) = 4H1 ( LQ, 

TI TL ( 16) =4HR )RH, 

TI TL I l 7 ) =4H0= , 

T I TL ( 1 3 ) = 4H l 00 * 

NTTL = 72 , 

XK=. 3306,-. 1734,. 10913, -. 03913, -.6 9° 98, .48287, 2. 6222, -1.40 68, 

£FNO 

£PERT 

TI TL ( l 3 ) =4HL , 
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XK=1.2583,-. 17087, l. 0046 , .0953 »-. 30686 , 1 . 0189, 5.9393 ,- U 734, 

SEND 

SPERT 

ICGNT=0 ? 

TITL( 18 ) =4H0 .04* 

XK=5 . 17, -.0 59 375, 4« 8719, 1 . 1309, . 106, 5.0654,7.2623,-2. 5309, 

SEND 


file PERT 2 INPUT 

IT CONTAINS THE INPUT FOR THE HALL AIRCRAFT 


SPERT 

N=4 , M =3 , IC0NT=1, IPRINT=1G00, IPL0T = 2, 

XLBL=4HREAL » 4H( L AM, 4HDA ) , 

YLBL=4HIMAG,4HINAR,4HY(LA,4HMDA) , 

NXC=l 2 , NYC= 16 , 

TITL=4HHALL,4H A IR»4HCRAF , 4HT CL,4H0SED,4H LQ0,4HP E I , 4HGENV , ^ H A LUC > 
4H PER , 4HTURB,4HATIU»4HNS P,4H=.l ,4H(MIN,4H K)T,4HAU= ,4H1.0 , 
NTTL=72,CNAMES=4H ,NCINDX=0, 

NS AMP=5 , 

P=.l, 

A=-3. 18, 1. ,-.06, .022,3*0. , . 0644 , . 63 , 0 . ,- . 27,-. 99 8, -10. 6 , 0. ,4.5, 

B=- 14. 4, 3*0. , 1.5,0., -2. 59,. 03 7, 2*0. ,-.96,0., 

XK =-3.97326316—02 ,-4.359845 1 E— 03 , 1.99 87654 E— 02, 1. 2020696 E-02, -9.9 E-OL, 

2. 5 59 335 7E+0 0, 2. 16 1 1822E-0 1 ,-3. 66 1 1261E-03 ,4 . 75 162 79E + 00 , -9. 4E-0 1 , 

-5. 4777920E— 01 ,-5.14106996-01, 

SEND 

SPERT 

T I TL ( 18 ) = 4H0 . 5 , 

XK=. 02469, -.001 394, -.Oil 66, .24 055, .0384, -.0438, -.30 16 5, .02299,4. 14655, 
.0114765, -.3113, 

SEND 

SPERT 

T I TL ( 18 ) = 4H0. 2 5* 

XK = 3. 5695368 E-0 1,-3. 9497081 E-02 ,-l .27873 79E-02 , 1.45003 13E+00, 2. 9E-01 , 

-a. 7618792 E-0 1, - 1. 29 22382E+00, 2. 6900148E+0 0,3.61001116 +0 0,-4. IE *-0 0, 

-1. 245 362 OE >01, 1 . 41 50043E + 0 1 , 

SEND 

SPERT 

TITL( 15)=4H(R0B,TITLI 16 J = 4HUST ) , T I TL ( 17 ) =4HT AU= , T I TL I 18) = 4HL.O , 

XK=-. 03489, . 1519, -2. 633, .137833, —.933, 2. 535, .2175, .002644, 5.5306, 

-.96 7 3 3, -.6965, -.525889, 

SEND 

SPERT 

T I f L ( L 3 ) =4H 0.5, 

XK=. 038 19 ,-.00 162 45,-. 012543, .4133, .03 86788, -.0443 5,-. 27592, 

.025 3 25 5, 4. 15 544, -.12295,. 007238, -.3 05 5, 

SEND 

SPERT 

T I TL ( 13) = 4H0 .25, 

XK=. 499 74,. 3 09 L3 l, — 3.344, 1 . 46229 ,-. 9 32 1 29 ,-5. 79063 , . 1 1 230 7 , - .Q ©// 37 2. 6 . 
5.979 13,-1.04167,-2.05667,-1.36633. 

SEND 

SPFRT 

T I TL I 15) = 4H{ LQR.TITLI 1 6 ) = 4ri ) , RH , T I TL ( 17)=4H0= , T I TL ( 13) = 4HIOO 
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XK=. 045 13, -005 74 ,.003798, .09 16, . 0222 , . 0 1 l 679 , o 0569 79 , .09242 , 
.03633, -.0754, —.02 08, —.009925, 

GEND 

GPERT 

Tl TL ( 18)=4H 1.0, 

XK=. 86 53, —.067577, .007833, -99586, -.058595, .016246, .1175, .93^ 89. 
. 35148, -. 55619, -. 27992, -.1, 

GEND 

GPERT 

ICONT = 0 , 

T I TL ( 18 ) =4H0 .04, 

XK=4. 8341 ,-. 42 64,. 02 78 7, 4.9 8, -.5015,. 01 23 7,. 4 1809, 4. 748,1.7639, 
-.73614, -3. L53,-. 86026, 

GEND 
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Appendix E The Sensitivity of Eigenvalues and Eigenvectors 


This appendix contains a more complete derivation of equations (16) and 

(17) presented in Chapter II. The derivation presented here was taken from 

Stewart (1973) , and more complete discussions of matrix perturbation theory 

7 

are available in Horscholder (1964) and Wilkenson (1965) . 

As previously described, we assume matrices A and A + E have eigen- 
solutions 

Av. = ,\.v. (El) 

and 


(A + E) vl = Xj vj (E2) 

with ( | E | l , = X] - and = v] - v^ all small perturbations of order 
0 < e << 1 . Using the notation and assumptions of Chapter II, we consider 
a specific eigenvalue .\ with corresponding right and left eigenvectors v and 
w satisfying 


v H v = 1 , w H v = 1. (E3) 

Given a unitary matrix [vU], transform A as 


1 — 1 
< 
cz 

1 1 

rc 

3>i 
i — i 

c 
cz 
1 » 

11 

v H Av 

v H AU 


U H Av 

u h au 


and v H Av = xv H v = x, 

U K Av - U H (av) = x(U H v) = 0, 


SO 

\ v H Au" 

0 u h au . 

The eigenvalues of U H AU will be those of A less \ since we have reduced A to 
block-triangular form by a similarity transformation. Next expand (E2) as 


[vl’.j H A [vU] = 
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with 

since 


(A + E) (v + Up) * (X + u) (v + Up) 
Ev + AUp * yV + .\UP 


EUP * 0( e 2 ) and yUp * 0 (e 2 ) 


,H 


Premultiply (E5) by U to obtain an expression for p; 

U H Ev + U H AUP a ,\p 

or p s (XI - U H AUf 1 U H Ev 

which yields equation (17) in Chapter II. 

U 

To derive equation (16), premultiply (E4) by v to obtain 

(v H A + v H E) (v + UP) = (X + y) 

X + v H AUP + v H Ev + v H EUP = X + y. 

Use the previous result (E6) to obtain 

u = v H Ev + v H AU (XI - U H AU) -1 U H Ev 


+ v H EU (XI - U n AU)~ 1 U n Ev 


- IlH/UlN"! ||Hr 


w = Cl v H AU(XI - U H AU) _1 ] 


.H' 


Ev 


+ v H EU(XI - U H AU) -1 U H Ev 
y = w H Ev + v H EU (XI - U^U)” 1 U H Ev 


v\ s l|w H l! ||Et| + ||EI| 2 j j (\I - u h au) _ 1 i ! . 


(E4) 

(E5) 


(E6) 


(E7) 


(E8) 

(16) 


To proceed from step (E7) to (E8) we used the representation for the left- 


eigenvector 


w H = [1 v H AU (xl - U H AU) -1 ] 


.H 


U 


H 


of A. 


(E9) 


, H- 


This identity can be derived as follows. Let [1 z ] be a left-eigenvector 
of the matrix 

[vU] H A[vU] , 


1 


u 

and solve fo.’ z . 


so 


H H 

Once we find z , y is given by 


y H * [1 z H ] 


,H 


[1 z H ] 


v H AU 

u h au 


= x[l z H ] 


[X (v H + Z H U H )(AU)] = [X Xz H ] 


or 

z H (xl - U H AU) = v H AU 
z H = v H AU(XI - U H AU)" ] 
which by ( El 0 ) yields (E9) and hence (E8). 


(E10) 
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